{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    z = 1 / (1 + np.exp(-x))\n",
    "    return z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialize_params(dims):\n",
    "    W = np.zeros((dims, 1))\n",
    "    b = 0\n",
    "    return W, b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 定义逻辑回归模型主体\n",
    "def logistic(X, y, W, b):\n",
    "    '''\n",
    "    输入：\n",
    "    X: 输入特征矩阵\n",
    "    y: 输出标签向量\n",
    "    W: 权值参数\n",
    "    b: 偏置参数\n",
    "    输出：\n",
    "    a: 逻辑回归模型输出\n",
    "    cost: 损失\n",
    "    dW: 权值梯度\n",
    "    db: 偏置梯度\n",
    "    '''\n",
    "    # 训练样本量\n",
    "    num_train = X.shape[0]\n",
    "    # 训练特征数\n",
    "    num_feature = X.shape[1]\n",
    "    # 逻辑回归模型输出\n",
    "    a = sigmoid(np.dot(X, W) + b)\n",
    "    # 交叉熵损失\n",
    "    cost = -1/num_train * np.sum(y*np.log(a) + (1-y)*np.log(1-a))\n",
    "    # 权值梯度\n",
    "    dW = np.dot(X.T, (a-y))/num_train\n",
    "    # 偏置梯度\n",
    "    db = np.sum(a-y)/num_train\n",
    "    # 压缩损失数组维度\n",
    "    cost = np.squeeze(cost) \n",
    "    return a, cost, dW, db"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 定义逻辑回归模型训练过程\n",
    "def logistic_train(X, y, learning_rate, epochs):\n",
    "    '''\n",
    "    输入：\n",
    "    X: 输入特征矩阵\n",
    "    y: 输出标签向量\n",
    "    learning_rate: 学习率\n",
    "    epochs: 训练轮数\n",
    "    输出：\n",
    "    cost_list: 损失列表\n",
    "    params: 模型参数\n",
    "    grads: 参数梯度\n",
    "    '''\n",
    "    # 初始化模型参数\n",
    "    W, b = initialize_params(X.shape[1])  \n",
    "    # 初始化损失列表\n",
    "    cost_list = []  \n",
    "    \n",
    "    # 迭代训练\n",
    "    for i in range(epochs):\n",
    "        # 计算当前次的模型计算结果、损失和参数梯度\n",
    "        a, cost, dW, db = logistic(X, y, W, b)    \n",
    "        # 参数更新\n",
    "        W = W -learning_rate * dW\n",
    "        b = b -learning_rate * db        \n",
    "        # 记录损失\n",
    "        if i % 100 == 0:\n",
    "            cost_list.append(cost)   \n",
    "        # 打印训练过程中的损失 \n",
    "        if i % 100 == 0:\n",
    "            print('epoch %d cost %f' % (i, cost)) \n",
    "               \n",
    "    # 保存参数\n",
    "    params = {            \n",
    "        'W': W,            \n",
    "        'b': b\n",
    "    }        \n",
    "\n",
    "    # 保存梯度\n",
    "    grads = {            \n",
    "        'dW': dW,            \n",
    "        'db': db\n",
    "    }                \n",
    "    return cost_list, params, grads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 定义预测函数\n",
    "def predict(X, params):\n",
    "    '''\n",
    "    输入：\n",
    "    X: 输入特征矩阵\n",
    "    params: 训练好的模型参数\n",
    "    输出：\n",
    "    y_prediction: 转换后的模型预测值\n",
    "    '''\n",
    "    # 模型预测值\n",
    "    y_prediction = sigmoid(np.dot(X, params['W']) + params['b'])\n",
    "    # 基于分类阈值对概率预测值进行类别转换\n",
    "    for i in range(len(y_prediction)):        \n",
    "        if y_prediction[i] > 0.5:\n",
    "            y_prediction[i] = 1\n",
    "        else:\n",
    "            y_prediction[i] = 0\n",
    "            \n",
    "    return y_prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztvXt8VPWd///8zEDIFRVBAyEkSBLCRQW0CCEXroEI2pt1bW2XtbTddguLK12tdqmtpRb9FSutdr91t9babWur1laBhCCEJBAuIkEFCeEaAyQIRCFXYDKf3x8zE3KZyZwzc85cP8/HI49HMnMun3My8/58zvvyegspJQqFQqGIHCzBHoBCoVAojEUZdoVCoYgwlGFXKBSKCEMZdoVCoYgwlGFXKBSKCEMZdoVCoYgwlGGPUoQQ9wshSk069ktCiFVmHNvNuf5FCLHNw3vpQggphBjg4f3HhBD/a+4IfUMIsVUI8Y1gj0MRnijDHsEIIXKFEFVCiAtCiCYhxHYhxGcApJR/lFIWhsAYg2bApJRPSinD3ngKIU4IIeYG4bxSCJER6PMqvON2JaMIf4QQg4F1wHeAvwIxQB5wKZjjigaEEAOklLZgj0MRvagVe+SSBSCl/LOUslNK2S6lLJVSvg99XRjO1de/CSEOCyGahRA/EUKMEULsEEJcFEL8VQgR427fbvv3Wb0JIa4TQqwTQpwVQnzi/H2k872f4phsnhNCtAghnnO+ni2E2OR8yjgkhLi32/GuF0K86RzTbmCMhnvxdSHEaSFEgxBiRbdj/UgI8X/O311um8VCiI+EEOeEED/otu1U57341Hmc51z3o9v1f1cIcRg4LIR4Xgixpte9eEsI8aC7AQoh5gkhapxPV88Bott7Y4QQW4QQ553j+qMQ4lrne38ARgFvOe/hw87XXxVCNDqPVyGEmODp5jj/n8ec//fjQoj7u733dSHEQef/bqMQIs35eoVzk/ec5/0nr/8FReCQUqqfCPwBBgPngd8DRcB1vd7/F2Bbt78l8KZzvwk4VvabgZuAa4APgcXu9u22f4bz95eAVc7frwe+CMQDScCrwN+77bcV+Ea3vxOAeuABHE+UU4BzwATn+6/geAJJACYCp3qPpdux0p3j+rNz+5uBs8Bc5/s/Av6v17b/A8QBtzrvwTjn+7cB05xjSgcOAg/2uv5NwBDn/lOB04DF+f5QoA240c04hwIXgXuAgcB/ADbXfQEygHnAIGAYUAE8223/E65r6vba1533exDwLLDPwz1KcJ57rPPv4d3u9eeAI8A453X/F1Dl7n+ufkLrR63YIxQp5UUgl6vG6qxzpXtjP7s9JaW8KKU8AOwHSqWUx6SUF4BiYLIP4zgvpXxdStkmpWwGfgoU9LPLIuCElPJ3UkqblHIv8DpwjxDCimOS+KGUslVKuR/HxOWNHzu3/wD4HfBlL9u2SynfA97DYeCRUr4rpdzpHNMJ4DduruNnUsom5/67gQvAHOd79wFbpZRn3JzzTuBDKeVrUsorOAxxo+tNKeURKeUmKeUlKeVZ4Bk35+6BlPJFKWWzlPISjgnsViHENR42twMThRBxUsoG5/8f4F+d13RQOlxLTwKTXKt2ReiiDHsE4/xC/ouUciSO1e0IHEbDE92NTrubvxP1jkEIES+E+I0Qok4IcRHHavNap5F2Rxpwh9Pl8akQ4lPgfiAZx2p1AI4VvYs6DcPovf2IfrZt7PZ7G85rFkJkOd1Ijc7reBLHStvTecAx6XzV+ftXgT94OOeI7vtKKWX3v4UQNwghXhFCnHKe+//cnJtu21uFEKuFEEed259wvtVnHyllK/BPwLeBBiHEeiFEtvPtNGBtt/9DEw4XUYqncytCA2XYowQpZQ0OF8lEAw7XisO1AoAQIrmfbVcAY4E7pJSDgXzXbq6h9dq+HiiXUl7b7SdRSvkdHG4UG5DabftRGsbbe/vTGvbpzX8DNUCm8zoe63YNLnpfy/8BnxVC3IrDnfF3D8du6D5GIYToNeafOY99i/PcX+117t7n/QrwWWAuDjdauuvQ7k4updwopZyHww1Tg+MJDxz/i3/t9b+Ik1JWebgORYigDHuE4gxArugWqEzF4YLYacDh3wMmCCEmCSFicTzqeyIJx2r/UyHEEODxXu+fweHHd7EOyBJCfE0IMdD58xkhxDgpZSfwN+BHzieB8cBiDeNd6dx+Ag7f/V80XWXf67gItDhXtN/xtoOU8iTwDo6V+utSynYPm67HcT+/IBw59/+O4wml+7lbcNzDFOA/e+3f+x4m4YgPnMcxAT/paYxCiBuFEHcLIRKc+7QAnc63/x/wqCvwKoS4RgjxpX7OqwgRlGGPXJqBO4BdQohWHAZ9P44VtF9IKWuBJ4C3gcOA2wIhJ8/iCCaec46hpNf7a3H4zz8RQvzS6YcvxOGTPo3DNfIUjiAgwFIc7pFGHE8gv9Mw5HIcQcDNwM+llL4UZn0Px0q4GceKVuvk8HscQVtPbhiklOeALwGrcRjjTGB7t01+jCOIfAHHJPC3Xof4GfBfTpfJ94CXcbicTuEIevc3mVtwfCZO43C1FAD/5hzXGzju/StOl85+HIF4Fz8Cfu88770oQgbhcOcpFAozEELk43DJpEsp7cEejyI6UCt2hcIkhBADgeXA/yqjrggkyrArFCYghBgHfIojINlfJpJCYTjKFaNQKBQRhlqxKxQKRYQRFBGwoUOHyvT09GCcWqFQKMKWd99995yUcpi37YJi2NPT09mzZ08wTq1QKBRhixBCS6W1csUoFApFpKEMu8JQ7HY7JSUlFM2ex+D4BKwWC4PjEyiaPY+SkhLsdpX1pwgtXJ/ZwsIiEhMHY7VYSUwcTGFhUdh+ZoOSFXP77bdL5YqJPGpra1k0bwH2T1rJax7CZIYSzwDasFHNOSqTmrAMSWBdaQlZWVnBHq5CQW1tLUVFi+hok6QlF5CaPIWYgfFcvtJGfeNe6hrLiY0XFBevC4nPrBDiXSnl7V63U4ZdYQS1tbXkTp3GouZkcu3JOHSseiKlZJulkXVJjWzbvTMkviiK6KW2tpbp03PJTrubMan5Hj+zR+srqKl7kx07tgX9M6sMuyJg2O12sm/KILc+ljw53Ov2FaKBqtRLHDx2GItFeQMVgcdut5OZmU3yNflkjOpX2h6AIx+Vc+ZiJbW1B4P6mdVq2NW3SuE3paWlyKZWcu39qfdeJc+eTGdTC5s2bXL7vvLTK8ymtLSUS22SMan53jcGxqTm095q9/iZDTXUil3hN0Wz55FcdpJ80V//ip5UcJrGmSMp3tLzi6L89IpAUDiviJZzKWSmeV+tuzhct5XEYacpLS02cWT9o1bsUUowVrvbd1Yx2XNDH7dMlkOp2rmjx2suP33uyVgea55AvhhBkojBKiwkiRjyxQgea55Abn0suVOnUVtba+RlKKKIqh3bSU2eomuf1OTb2FEVHj1GglKgpDCH3qvdVUxxrHbbbVSXnWTpnsWmrHZbO9qJ1/lRimMArZeu9p2w2+0sKlzAouZkh5/eba8fEEKQJ4cjm+GuwiLlp1f4RHtbKzED471v2I2YgXG0tbeaNCJjUd+ICCGYq92E2DjasOnapx0bCYPiuv422k+vUPRHXHwCl6+06drn8pV24uMSTBqRsSjDHgH0Xu26S9uCq6vdhc3J3FVYZJhbZsa0HKo5p2ufanGOnGnTu/5e+/QacpuHeBx7b4QQ5LUM4dmnfq7rvAoFQE7ODOob9+rap77xXXJyZpg0ImNRhj0CCOZq1263M2N2AX+3nuC7soIlcgvflRX8Ur7PB/I8djfBeSklFQlNPPjI97peM8pPr1Bo4aGHllPXWI7W5BEpJScaynloxXKTR2YMysceAfiz2p0/f77P5+3y6Te18rnO9D4ZLK9zlFc4zDJ5C8niqj+z0tLIgOsTmTdvXtdrRvjpFQqtFBYWEhsvOFpfoSmP/Wh9BfGJ1h6f2VBGrdgjgGCsdnv49Fvc+/Qf5zPMZxRPsZdG2eZYqYsG1ic18lZpcY+gpxF+eoVCKxaLheLiddTUvcmRjzyv3KWUHPmonJq6N9mw4a2wCdSHxygV/RLo1a4en36+GMHnGM3PqeanifupSr3kVk7ACD99uKEKsYJLVlYWO3Zso/FCBZt3/4TDdVvpuNSM3W6j41Izh+u2snn3TzhzsTIk5AT04LcrRggRC1QAg5zHe01K+bi/x1VoJyE2jrZ2G0nEaN7Hn9XuVZ9+use0xO7kM4JN1lN8/fvLefTRR92uepY/vIKlexaT1+x5ouiOlJKKxCZ+/Uh4thMNVmqqoidZWVkcPlzDpk2bWLPmWTZUvkZbeyvxcQlMz8nhxdXPMW/evLBZqbvwu/JUOL6FCVLKFmdX9m3AcinlTk/7qMpTYzGy8jNY54smvRklmKbwlYBVnkoHLc4/Bzp/VIfsALL84RVUJjXpivD3zkrRgxk+fYvFwrrSEtYlNVIhGvr1eXry04cDwU5NVUQHhnwrhBBWIcQ+4GNgk5Ryl5ttviWE2COE2HP27FkjTqtwUlhYiGVIAtssjZq2d5eVogezfPpZWVls272T7akdPJl0gAp5mmZ5GZu00ywvU8Fpnkw64NFPHw6oQixFIDDEsEspO6WUk4CRwFQhxEQ327wgpbxdSnn7sGFee7EqdBDo1a6vGSzxMbFeg4VZWVnUHDvC86++TOOskayMq+Y7lgpWxlXTOHMkz7/6MgePHQ5Low6qEEsRGAxXdxRCPA60Sik9fhKVj90camtrWVToyCt3KSPGMYB2bFSLc1QmNmEdkshbpcV+GUZffOxvyRNsGdTI0JjEqFZtHByfwKr2KSQJ7YHuZnmZlXHVXGhr8b6xIqIJmI9dCDFMCHGt8/c4YC5Q4+9xFfoJ1GpXr0+/QbaygTo+fznNVB2bcEgfVIVYikBgRFbMLcDvASuOieKvUson+ttHrdjDGz0ZLHYpeZgq7iKdApHi9di+ZruEi467WrEr/CGQWTHvSyknSylvkVJO9GbUFeGPHp/+K+IIMVjJR5vbxpdgYTjpuEdjIZYi8IRXrpgiZNCawfJuTBMLGGVasDDc0gcDnZqqiE6UYVf4jBafvk1IU3Vswi19MNCpqYroRKk7KvzCYrEwf/58jyqRrZfMDRYGS9nSV1xurNyp05DNjonGU+VppaWR9UmNbCvdGXaFWIrgoj4tClMxW7UxHHXco6EQSxFclGFXmIrZwcJwTR+M9EIsRXBRrhiFqZit2hhoZUsj8ebGUih8Ra3YFaZidrBQpQ8qFH1Rhl1hKmbr2Kj0QYWiL8qwK0zHzGChSh9UKPqiDLsiIJgVLIwWHXeFQg+GqztqQWnFKIwmUMqWCkUwCZhWjEIRCqj0QYXiKsqwhyjhIEEbarjSB4u3bOJCWwu2zk4utLVQvGUT8+fPV+4XRdSgPukhSG1tLdmjM1h672KSy06yqn0Kv5EFrGqfQnLZSZbeu5jsmzKCqlJoBGryUijMQfnYQ4xo6WAfLvrpCkUoodXHrgx7CKGngQX43pQi2ETL5KVQGI0KnoYh4SZB6wvhpp+uUIQjyrAHEG8+5Wef+nnEd7CPhskrVHB93goLi0hMHIzVYiUxcTCFhUUqhhHhKFdMgNDiUz7VfJ4VTGKMuEbzccOtH2bR7Hkkl50kX2hrlQdQwWkaZ46keIsy7lqpra2lqGgRHW2StOQCUpOnEDMwnstX2qhv3EtdYzmx8YLi4nXKzRVGKB97CKHVp1zOaf7BcR5hCskiXtOxbdLOdywV2Do7jR62KQyOT+CJ9sl8RAtlnOIQn9KBjVgGMJZrmUUKExiCpds9CoXJy263U1payjPPrKWqajvtba3ExSeQkzODhx5a7pA2CJE4R21tLdOn55KddjdjUvM9ft6O1ldQU/cmO3ZsU8Y9TFCGPUTQGxAtl6fYxEmeYGoP4+aJUDB6erAIQTLxDMTCbEb2eXLZwkmuYGcZt3RNbsGevLSsfhEdjBmVyr7336O1o52E2DhmTMth+cMrAmr07XY7mZnZJF+TT8aoAq/bH/monDMXK6mtPRgyE5PCM1oNu9JjN5mrPuV00OA6z2cEZZziQ5qYyPVetw8nCdra2loGYmU+qeQxosdKMokY8hlBnhxOJQ08xV4ekY4nl2Dqp/dY/U7sufqNHZREZloBGaPyOVy3lT3v/JmH7BNIJ4m2dhvVZSdZumdxQNM2S0tLudQmGTMxX9P2Y1LzObG7nE2bNild+AhCTdEm40tPzlmksIVTXrcNJwlaVzbMl8kkX6T0mw2TL0bweW7iOT7ALmXQJi+73U5R0SKy0+4mY1RBv2POSp/FpIlf5n+tRxAIkkQM+WIEjzVPILc+ltyp0wJSUPbMmrWMSvY81t4IIUhLzmfNGm2NTRThgTLsJuNLT84pDKOWT71uF04StK4nl3y8u6MA8hjOAAQHOB+0yatr9ZuqbfWbkT4Te2wCH9LU9Vqg0zardmwnNXmKrn1Sk29jR1WVSSNSBANl2E3G156c7dgiSoLWlyeX2YzkbxwP2uTly+o3I2M+pdYzfd4LVNpme1srMQO1Bd5dxAyMo6291aQRKYJB6FuEMCchNo42bADYpeQDeZ5fyvf5rqxgidzCd2UFv5Tv84E8j91pyNuxMUBYI6qDvS9PLpMZymlagzZ5+bT6HX47R+QnfV4PVM1BXHwCl6+06drn8pV24uMSTBqRIhgow24yrp6cjbKNlezidY4yiaGsZhovMJPVTGMSQ3mdo6xkF42yjWpxjjkzZ0WUBK2vTy6dArfXGYjiG19Xv5fsl92+N1kOpWrnDr/H1R85OTOob9yra5/6xnfJyZlh0ogUwUBlxZjM8odX8K3d9/NG6zE+z03kMdxrNkhSXCL/+8izEdXBPiE2jrZ2G0nEaN6nHRuJsX2zYXqnHy7MW90j/XDJA8sMKb5xrX5jByVp3ufylXYGWWLAjRctjgG0Xmr3eTxaeOih5Sx5YBkZo9znr/dGSsmJhnJ+t/p5U8elCCxqxW4yc+fO5XxHC59nNPlihNdskM8xmk8utTJnzpwAj9RcXE8uenCXDeNKP0y+Jp85U1eSmVZA7KAkLBZrV/rhnKkrSb4mn+nTc/3KRPFp9duwhwxxndv3ApG2WVhYSGy84Gh9habtj9ZXEJ9o7RHDUFIE4Y8y7Cbz9ttvkxw7mDy0ldDnM4IbYpPYvHmzySMLLMsfXkFlUpPHgHBv3KVy6kk/zBhVQHba3dx5510+G6KHHlpOXWO5rjEfPlJCYeeNbt8PRNqmxWKhuHgdNXVvcuQjz2OXUnLko3Jq6t5kw4a3umIYtbW1ZGZms+SBZbScTWFh3mq+sui3LMxbTcvZFJY8sIzMzOyw7wUQ6SjDbjJrn15DfutQXZkVBW1Dw0rYSwuFhYVYhiSwzdKoaXt3qZx60w/HpObT3mr3ORNF7+r3yImtWDvaGM+QPu8FsuYgKyuLHTu20Xihgs27f8Lhuq10XGrGbrfRcamZw3Vb2bz7J5y5WNlDTiCQT0MKc1GG3WR8ygYJQJBNK0Z1ObJYLKwrLWFdUiMVosGnVM5AFd+4rnnBgoU0NJxm13t/4NDxLf2O+fDxMj7Y/wrLO8e5lYIIdM1BVlYWhw/X8OJLz5E47DQbKh/lT+u/yYbKR0kcdpoXX3qO2tqDXUY90E9DCnPxO3gqhEgFXgaSATvwgpRyrb/HjRR8zQYxO8imhd6KlKuY4tB18bFcPisri227d7KocAGVTQfIax7CrVzPcS6yyXqGI/ITrkgbVjGA6RNzOHbsGBkZGV3GvWrHdhbmrdZ1DanJt7Gh8jVd19w9MHtXwdN0XLrI1t1rOXT8bcaNKSQ1+TZiBsZx+Uo7HzXsofZICQM7WnnMfmsf8TYpJZWWRtYnNbKtdGdA0zZdPWC1BOCVFEFk4bcImBBiODBcSrlXCJEEvAt8Tkr5oad9okkEbHB8Aqvap5AktGeDhIKwl5ldjux2h3vkp4//mJ3vVhM76BomZi7yKi1rtVj5yqLfYrFYNV+H3W7jT+u/SWenTdM1e1JFlNLO6bP7OXTsbRrP12CzXSI+PoEpk2/j2OEaYtvs5Ldcz2SGdhWYVYtzVCY2YR2SyFulxSGdnlo4r4iWcylkpnkXDnNxuG4ricNOU1pabOLIFN0JWAclKWWDlHKv8/dm4CCQ4u9xIwWjskECidldjiwWC6NHj+bA4SN85ub7+dycpzX5c80svvHmihDCQsoNtzB72kN8ZeEL5Ez+OsOHj6C8Ygv1Daf59Wt/CLmaAz3ZLUqKILIwNI9dCJEOTAZ2GXnccGb5wytYumcxec2eDWR3pJRUJDbx60eCJ8qkV5Eyz55MZdMBzY/lvY2oJ1z+XIA777yLnOmO9EM9q0qtxTf+uiJCreZAb66/kiKILAxz+AkhEoHXgQellBfdvP8tIcQeIcSes2fPGnXakMeIbJBA44uui55yeV+zWwpm5upOPzzRUM5DK5Z73TaSVBF9yW5RUgSRhSGGXQgxEIdR/6OU8m/utpFSviClvF1KefuwYcOMOG1YYEQ2SKAxO5PHVyO6dWul38U3nogUV4Sv2S3Tp+coKYIIwm/rIRyfnN8CB6WUz/g/pMjDlQ2yPbUjLIS9zM7k8dWI7tyxw6/im/6IFFeEr09DM2fmmfY0pAg8RiwLZwBfA2YLIfY5f+404LgRRVZWFjXHjoSFsFd3RUqt6CmX98eI+lp8441IcUX4+jRUvnWbaU9DisDjd/BUSrkNTSE2hZ684mAyY1oO1WUnydcogwD6Mnl8FddyGVFX8c2mTZtYs+ZZNlS+Rlt7K/FxCUzPyeHF1c8xb948Xe4sly6MlsCslHZOf7yf6oOv0XGpA6vFGjKNrf3J9X937ztMn54LoLkJdjj0AohGlLqjog9mZ/LoMaIuevtzjZ4ktaoiXmxpoGzXWqzWgYwdPbdP7r0nZUm73U5paSnPPLOWqqrttLe1mjIZGPE0VFS0iBMN5aQl5/coxqpvfJe6xgriEiy6noYUgUcZdkUfujJ5WhvJk95b2enN5AlFadnuujCeUjAvtjSwcduTTBp3T5+xd29sfbS+gunTc7uMX6BkhiE0n4YUgcfvylNfiKbK03DFVXm6sDmZvH4qT7vK5XVWnmZmZpN8TX6/eewujnxUzpmLldTWHjTVoHirPH1zy2OMG7OArPSZmse8bt0/mDEj3+0xrx67p3vDH+NeWFhEy1n9FaRJwxrYWLrB5/MqAkPAKk8VkYmZmTz+SsuaRX+B2ROndiOERbPBHJOaT1tLJ3NmFwZUWMsXqWGV3RJ5KMOu8IiZmTxmZbf4iydVxN0f/J5xY+bryja5JmEM7W32gMkMgzGNNhThj3LFKIKKSxBszZpn2VFV1cOfu2LFgyHjz01MHMzCvNW6fNdv7/j/SBsxNeDCWv25lFwY6f5RBA6trhgVPFUElXBJAfUl2+Rs0xFyp3xb1z56ZYbdobJbFMFfCikUYYAvBUw2W0fQqln1NtoINVTfVf9QK3aFQgO+5N5brTF+pR76S7g8DfUmkOmhkYpasSsUGvAl2yQmJk4JaznRugI3s+9qND0FKMOu0IVRPVDDDV+yTa65NsH01MNwMFa1tbVkZmaz5IFltJxNYWHear6y6LcszFtNy9kUljywjMzMbGpqakzru6p1DJHSoFtlxSg007sH6mSGOnqgYqOac1QmNenqgRpu6M022b69goUL7zatEKu3y6K3vMHBY8W0tjU52vglJAZFy0bPPXv/8GskxV/P3GmPa65I3rz7J7z40nP9upsiKUtIa1aMMuwKTZjZAzWc6GlMPWebbNjwVpecgBlGRetxD9eVs6/mdWZNXc6nzaf69JE1E70VxuvLf0hW+hxD00NDtcrZV1TlqcIwjOyBGu6uHL3ZJmYUYulpppGVPpPJ475I1b7fkjEqz2f/tC/o1Ya/2PKx4c1OfNWn96dILBRQWTEKrxjVA7W3K2cVUxyunHYb1WUnWbpncVi4cvRmmxgtrKW3P2vGqAIOHd/M6bP7Sbnhlh59ZM1cmerVhrfZ2g1PD/Wn5WG4ZRN1R63YFV4xogeqy5WTezKWx5onkC9GkCRisAoLSSKGfDGCx5onkFsfS+7UaRETxHLhmgxKS4tpbrlAZ6eN5pYLlJYWM3/+fF3G1RdjNXb0HA4d39z1WiBWpno7ZQ0YEGd4s5NIaXmoF2XYFV7xtweqka4che/G6uPzh7r+DkQzbr3VujcOHWt4emiktDzUi3LFKLzibw9UPa4cu5RcKwfy8cnjJMUl0HHlEgmxccyYlsPyh1cEtTtRqODOWLm6Oh06sZkz5w5hs7UzYEAcNw4dy9j0OSQPzeaKraPHPkbIF/SHXm34selz2PvhXw3V6fdXnz5cUYY9xHF13ln79Bq276yitaM94IYuITaOtnYbScRo3qd7D1StrpxG2caveJ+BWPhcZzqTO4eGpQ/eDLp3YEJYehir3l2dciZ9o0faY/XBV7HZLjHA2vP/Z/bKVG+17ogbJrJj34scrttKVvosr9trUaY0oltXOBLdS58Qp7a2luzRGSy9dzHJZSdZ1T6F38gCVrVPIbnsJEvvXUz2TRmm+6NnTMuhmnO69uneA1WLK6dRtvEUe5nPKB7nM3754MOhaEcPvYtrkoeO63JZuLo6jc8oYmHBE24rNRcWPMH4jDu7tndh9spUb7UuCBISEvjw+N8N0+mPVn16tWIPUXrmjaf3WO0mEUM+I8hrHs621kZyp04zNW/c3x6o3lw5din5Fe/zeW4iX3huoO3ywctmuKuwiIPHDvf5UkeazkiPfPWJDhdFfOx1VB98lTGpuZTtWsukcff0uyJ1pT2CpGz3L7l71k8RwmL6ylRLu8HuHK2vYPC1caxbV8rChXcbokzpyxgiQZ9erdhDkFALNnb1QLU0atq+dw/UhNg42rB53P4ATcRgJQ/v/VXBkU7Z2dTSJ6PDTJ2RYOApX33EDRPp7LzCngOvYLXGkDFKW9pjZtpMrJYBnD6739CVqacnpAULFvLoow9z8MQ/dK3As7OzDVOmDNVuXWYTFZWnoeCn1kNJSQnL7l3Mo80TNK+Qn0yh5hgaAAAgAElEQVQ6wPOvvmxa7q2nHqh2KTlAE2Wc4hCf0oENKxZy7pjGoz9aSWFhIQvnzie57KTH1fgv5ftMYmi/q/XeVHCaxpkjKd7iMO6RVmEIjs/BNx5YxuypK/t8Di62NLCh4glum3Cf7krN+sZqRg2fYsj1e5M1qGssZ8DATuzSjt02UFO1rhnorRgOVVSjDSfhWBTjT9640Ya9+6TY1tHOn+yHKOYERTKNFBJ4kYMMwMIcRvIA2Ve1Y3adY+m9jnv78H89xuo9P/DoyjnEpzxAtq5xTZZDWelMpwT9RTtjUvM5sbu8TxFVKNFfvvrgxOFIafcp7fGd/X+i7fJJduzY5rdR7+0mcuF6QsoYld8ll/D00z/jL395ze8CLV8wukgs1InoFXug9E2MfiIYHJ/AqvYpJAntWSjN8jIr46q50Naie/yecCf6FYuVas7xJic4Rzv3kUkBI7ze2/iEBGY3JpEn+7pblsgtvMBMrEL7PbJJO9+xVGDr7ASgcF4RLedSAt6Gzky8teP7wz8Wc/9dL2KxWDUf02638X9vLeHQoRq/FjKR+IQUDkS9Vkyg/NRmZK74mzduBJ4qRQcKK7dzAxLJV8hkpkjRdG8tCN5KbKBCNPTxc8Y6V/l66J5OCZFZYeituMbXSs3EhCS/n06jVYMlXIhYw361KCZZ0/aeAnL9YVaZvLdgozt6Gzp/8DYpXg12avOJ59mTsTRfYvUzP2d7agdPJh2gQp6mWV7GJu2MYbBf6ZQQmRWG3trxmVGpqRV/NFgU5hOxht0IfZP+MPOJwN+8cX8VFL1NimWcYhaeV+q9cd3b1/74CjXHjvD8qy/TOGskK+Oq+Y6lgmMxrZRaT+rKNa5IaOLBR77X9ZovPUlDvcLQVVzjibHpczh0/O2g5GhH4hNSJBGxht1ffRNvmPlEsPzhFVQmNflk6IxwDXmbFA/xqc/31iWGVbxlExfaWrB1dtLc3kbCyGE+p1OCdyPojlCvMPRWXONKezzykfauTkblaEfiE1IkEbGG3Ww/tZlPBL7mjaelpRniGvI2KXZgM/TeWiwW1pWWsC6p0a0P3oWUkgrRwPqkRt4qLe4RhIvECkNv7fiEsDDrjuXsO/gah+u26s7R9qdCNxKfkCKJiE139FffxBvbd1axCn2Por1T9DzhMnS5U6chm+mRN94dKSWVlkbWJzVSUVLF3Qvu7HINeRLb0lK96W1SdAU7jby3WVlZbNu9k0WFC6hsOtCVhRPHANqxUS0crfesQxLZVto3eykSKwxdxTXTp+cCuO2UNDhxOIUzHmPTjqc4cGQDEzLu1FSp6W+FbrRqsIQLEbti99dP7Q2znwhchs5dsLFZXqaC0zyZdICq1Ets272TEydOGOYa8ha8Hcu1ptzbrKwstz74lXHVNM4cyfOvvszBY4fdGppIrTDU0oHpnQ//h+Th1/PMs6s0VWoaUaEbiU9IkUTE5rGXlJSw9N7FPKajevOniQf49WvaqjcDlWtutztSxJ596udU7dxB66V2EgbFkTNtOg8+8r2uooqi2fP6re50R+/qTRfejvWBPM/rHOVxPmPKvfWHSKkw7I3rc7BmzbPsqKrqUVyzYsWDmotrjMo/V3nswSGgzayFEC8Ci4CPpZQTvW0fCMNut9vJvimD3PpYt0UxvakQDVSlXnLrmnCHkYbUCIyYaFyFVo9//wdUv7cPG3ZiGcBYrmUWKUxgCBYhsEvJSnYxn1Garl/vvfUXo4xgJNKfTIE7pJRs3v0TXnzpuT6TslmNuhWeCbSkwEvAc8DLBh3Pb3zxU28r3an5C++v4qHR+Osa6l1lej85V+UBOMfrHOUVDrNM3kKyiGeZvIWn2AsS8nB/D3y9t/6itydpNOFL/vmoG/NYtPCzDIqNISdnBg89tJzCwsIuN1FR0SJDlBgVxmHIN01KWQE0GXEsI9Hrp9bzwfNX8dBo/Clq0lJo9TifYT6jeIq9NMo2kkU8DzOZNznOY+zy6976m3ev0I4v+eejht+OxWJlYd5qWs6msOSBZWRmZlNbW9ulwWKEEqPCOAzzsQsh0oF1nlwxQohvAd8CGDVq1G11dXWGnFcLWv3UevGkeNidHqtWEzXTfXUNNRSkcOz4Mc0uq3J5io3UM59UtiV9guW6BB754Q947Y+v+HRv3enRdH9SqExqCjmRtnDGarHylUW/1a0v88d13+Brd78EKPdKMAk5dUcp5QvAC+DwsQfqvGDeo7m/KXpG4qtr6Otzv8rRpz7Q1I8UIJ8RlPAR+2+J5/mn13YZ7iVLlugecyg1E4kWfO0BOnBAbNffQoiugOmdd96lAqIhiPpv+ElGRgbP/vo5rh2dwl8tR1nBdv6VrTxs2cmHt8Tzq7+85DFFz0h8dQ1te7tMd6FVkUhj2PXXM3/+fJ+/0KHWTCRa8LVC94brx/Z5XQl7hS7KsPuBq3z/3+97gJvfb+dn9jv4bwp4lly+bM/g/LGTLP/Odzly5IjpY/G1enP7rh2mSi94IhAibYq++JJ/XnP8bcaOntPnPSXsFboYYtiFEH8GdgBjhRAnhRD6n8vDDLOUHXujJ7DoS7A4WBLBZou0KdzjTaagN0c+Ksdu72TEMPdZzErYKzSJ2AIlMzE7R96FL4FF10Tww0f/i4P793PZbqMTSYxlALfefAs/+tmqHi6UYDX1CJVmItGI1vzzIx+Vs+/g68zPfYzBie4/53a7jT+t/yadnfoyshS+EXLB00jiqhshXVPAMc+eTGXTAV1t2HwJLAI9JoJ/5o6rE4H9HJVO11D3iWDGtByqy06Sr1FbHfRJL3giFJqJRCta8s8PHd9Mp93Wr1EHJewVqijD7gN63QgSGN0cwwP3fZWWS21eW+f1DixqEfSaP2sOra2tujNMglVoZbZIm6J/3PUAbW1rxmodRPLQcUwe/yVGDJuI8NKuUAl7hSYqeOoDerTeG2UbK9nFIT5l/qdDNemj6w0szrDfyMcNZ1h4UX+GSbAKrcwWaVN4x5UGXFpaTHPLBTZs2MANQ0cya+qDpNxwi1ejbrd3UnOimPPnz5OYOBiLxcrAgYMYNCgBi7CQkJCkSQJYYTzKsPuAVjdCo2zjKfYyn1H8iM9oDrDqfSL4kE8YImPIk/ozTIzQQvcFf5qJ+KMjrvCMnsDqxZYG/vb2CuydEGO/mYV5q7l/0W/5wtxfcNv4L3PdNaMYIBI5W5/Yo1JVERhU8NQHtAT+/BHKujYxSVdg8ZfyfSYx1C9BstraWhYVLsDe1Oq+0CrRUWj1VmmxITn5vgag/1GynoUL7+7SEU9NntJDR7yusbxfHXFF/2gJrF5obqC48gmmjL+XzLSZ/QRfK9h38DUKZzzG2U9qVaWqAWgNnqoVuw9ocSNcbfjs3WhBz1W03sCiP63qXPijhe4Lvjwp/PI3v2bGjHy/dMQV/eNN/732RBkl2xxGPSt9Vr9uv8y0AiaN+yJb3/klY1LzyE67mzvvvEs9UQUAFTz1AS0BR18bPj/71M91BxaNalUXaFVEPZIMFSVVLFx4N9lpd/er/63K3f3HXWDVJX2cmZXJtdfcQGbaTE3HyhhVwKHjmzl9dj9jUvM5sbtcV3aYwjfUJ94HtAQc/VlF6w0sulrV6SFUMky0PimcOHGCS22SMan5mo6ryt39o3dgtbPTRnPLBYZeP4ybUmbrWrCMHT2HQ8c3q0rVAKIMuw9ocSP4s4rWG1jM4pqwzjBxGZHiLZu40NaCrbOTC20tFG/Z1FVM5YuOeFpyPmt+bowRUdLCDnyR/U1Nvo2Pzx/q+l1VqpqPMuw+4q18Pwarz6tovSmIw0Q8JaLepwwTvQTLwPlqUMrKNvvta3dpAi29dzHJZSc1paxGKu1trcQMjNe1T8zAOK7YOrp+b2tvNWNoim4ow+4H/bkRhlx7nc+raL2BxX1JzcQOH2J6LrpWA1dTU2O48ffVoHR22vzS6gmUJlC44JL91UN32V9VqRoYlGH3E09uhBf//Aef87RBn6DX9nd2sbHsbVNz0bUauBn1sUwafzPfuud+Q1e3vhqUQZYYnyV/w1Fa2Owcf39lf1WlamCIKsMeSDeCERWdelIQzWwDqMfA5cvh3CfH0NnaQS7DDVvd+mRQGvaQIa7zWfI33KSFa2tryczM5uv/spT6o5JrEtKxWgfR1tpC2ZYyvviFL5OSkkpNTY3P5/BH9ldKyYmGch5asdzn8yu0ETUFSsFowRaM1nlmtAEsKSlh2b2LebR5gmY9mR/zDvcwhonierfb6FW8LCkpYckDy5gzdaXmMRRv/j7/3HoDE8X1fQqytOBru0G95zECV2FR2o0FHKuvwmodyNjRc/sUcB08upG2jvO8+dYbPqUc2u12MjOzSb4mv9+0UxeH67by4dGN3D3rpxytr+TMxUq3Kah2u53S0lKeeWYtVVXbaW9rJS4+oUfzbJW2qr1AKSoMe0+lRM8GdpulkXUG9yYNdEWnGfhk4ORp9nGOfxe3uH1fSsmTSQd4/tWXNRkY3QbleBnHDvyNn3ZOwSKET5K/4SIt7Lo31wy6lZpjm5g07h4yRnmW4z1cV867B/7E+x9Uk52drft8emV/vVWe1tbWUlS0SFUTa0BVnjoJtp/UrIrOQLqV9IieuZjMUGr51OP7ehtnWCwWiovXUVP3JoeOb+k3jnD4eBkf7H+F5Z3jsDj/375I/oaLtHBpaSkdrXaO1Vcxadw9ZKZ5TgsVQpCVPpMpE77MnDmFPn1OtFSnri//IQeOFDN29Fze+fB/OHOx0qNRnz49V1UTG0zEV54GQjvdG0ZXdPZ2K61iisOt1G6juuwkS/csNtSt5KuB66Cz320my6H8144qSkpKWPv0GrbvrOpX0thlUCZOuIWjR0vJzFzQU0e8YQ9HjmzE0tHKY/ZbSRZXs2h8KcgKF2nhZ9asZXBCBh0dh8gYpa2AKyttJoeOb/T5c+6+OrUFq2UgFouVy1c6iI9LYFhqK6uffc6tC9But1NUtEhzNXFnZycTJ9xKwcwCVqx4ULln+iHkXTEu35uWL747wslPqgWtbqUKGvgzh7lCJ4lx8Zrvlzt8dUk8yk6eE54NzSnZwir2MDJpqK64x4JZc7FvPUSD1cYR+QmX7JcZZIkhQ1xHYeeNjGdI10rdRST72BMTB3NNQjppI6aSmebdTeWi9kQZSTc0UFpabOLoPFNSUsLXF3+XudMe1xw3Wbd1JclDx3Gh7WhUumciwsduRMAzXPykWtCriFguT1FKPQ8zmfc473OA2Awfe6NsYzV7+RyjKWCErrhHSUkJS+9dzGM6grk/TTzAr1/T5s93Eajz+IvVYsVqHcTn5/6c2EFJmvfruNTMhspHaW65YOLoPJObOxN7a4auyehw3VbqG6uZNfVBjtZXRJ1iZNj72I0qDDHSTxrssnK96Xf5jGAgFupp8SvV0Bft9C2cZDYpbt+3S8mveJ8vMJqZwrNQmqe4R6CagwSrCYle4uITsNk6fCrgClYVqN1uZ+fOKp/lCVzuGaUY6Z6QNOxGBjwTYuMMEcgKhbJyvQ04hBDMZiRbONX1ty8BYt0GjgZsSMYzxO37DkljC3ka+6zm2ZNpOfkxiXHxWC0Wrk1MIjk5mTdiT5raHCRYTUj0kpMzA6s1xqcCrmBVgZaWltLZafNLngCU2JsnQtKwG1kYYkQLtlApKzcqO0VvIY0eA1cuT/EGx1jKzX383C4cksYj3U5Qdin5QJ7nGcuH/JvYxhJZxnct27EwkBsuD+S/ZT6r2qcwZtcnJImB/FUcZVXCB5oLsvRWZppZ+GUUDz20nJiYOOob9yKlnVNn3mfLrl/w5/Xf5g//WMyf13+bLbt+wakz7yPl1esLZhXoM2vWMsDHycglT+BAMjj+Jr72tQdUN61uhKSP3ciglb9+Ul87/WgtvNGD1WLhN7IAq5delN2xSTvfppz/FbN6jtOHIJ+3nPzy+HM0tn3KAkaxSKZ59Jv/GxU8zfQ+cY9G2cZa64fYYxPJyCzqk8984PA64jraWd45nmQR7yjuEg28EXeSW2++hX3vv9dvQZY/+dJmFH4Zhd1uJyUlldbmTizC6rE46dDxt+nsvMKsO5aTlJDM27ue4He/fz4o2uj+BHxPntnH7Dv+g4stDZTtWovFMoDsm+ZFRf57WAdPjQx4+muYfam61FN4owcjs1N8DRB7M3BpaWncveDOfguyjjV/zAvM7DFBNco2nrS8x80330dGf+3W6rbywQev9Ehn1DKZai2qCdeAXHFxMXfd9Tk+M/GrHvPYu7erG3vTXJovfxC0RiRWi5WZUx9kX83rLCx4QvN36++b/5Opt/wzSfHD2LjtSa/FWOH6//REWAdPjQx4+usn9cWvrafwRg8+uZU4RxbX9nnd10Iab9rp2dnZXguyEmPje8Q97FKy1vohN998H5ne2q2lz+Lmifex1noQu/N/6c211Dtfur/jh2NAzm63s3Tpcu64+WtkpbufFOFqu7pbs7/Ah0c2sG7dP4L2lBEXn8D1146ms/MKRz7y3jwbHBkx7R0XGT50PGW71moqxgrH/6cRhKRhNyrg6cIfP6lPfu1e/USNwsjsFDMLabwZ/xnTe05QB2jCHptEhoZ2a1LaiYsfwkWrnX+zXPXBX2q7wg8e/r7bL29paWlEd1/quj6NxUmZaTMZct1w6urqTB6ZZ3JyZnDyzD5m3bGcfQdf43Dd1v6rieu2svfDv3L9daNpOPshVmuM5mKscPt/GkFIGnYjAp698bW0P5TKyo3MTglmB6XeE9Qm6xkyMhd4fSq62NLAm1seY9/B15g0/l4+V/gs99/1IncX/oLUm7/A4aMNZGZm9wle+9x9KUxauPlyfaNHzAzq9blUIpMSkpmf+xgfHilhffkP+8gTHK7byvryH/Lh0RJiBw1mYuZCDp3YzNjRcyL2/2kEIWnYfVmZ9u4I5C7n/NrEJJ596ucsf3gFn7RcdNuCrTdGPz34gy63kjztMTvFnw5KRtB7gjoim7zmM19saWDjticZn1HEwoIn3GqKfHbO0241RXztvhQuLdzC8foKCwuJjRccra9gcOJw7p79JJPHf4n6xmr+vvk/+eO6b/D3zf9JfWM1k8d/iXE3zUcICyOGTeTMuUNhd72BJiS1Yrq++K2NmgKevQtDjNRSmTEth+qyk+RrzLkGc1fDLrfSosIFVDYd6Bug5BxbOIkNySNM6aGX4iJYhTQuXBNU7tRpyGa4ZL/Sbz6zlPYePlVPuHyqAHfeeVdXYNBd9yUp7Zz+eD+HTmzmzLlD2GztDBgQx41DxzI2fQ7JQ7PDpoWbr92lgnl9LlG36dNzAYe7JOWGW0i5oWelsksl8r2avzE/9zGEsGCztYfd9QaakDTsvb/4mrTMS3disVh6aamk99gviRjyGUFe83C2tTaSO3Wa19zj5Q+vYOmexeQ1ey6U6j2misQmfv2IeY99LreSKztl5c4dtHS0Y5UwggS+yE1M4Hq3K/Xe9ytYdJ+gLPUDuHylzWM5/OmP9+v2qZ7YXd4lcOXqvuQ6vitNzpUWmDPpGz3S5KoPvorNdonYQbFezhQa9L4+LQS6OMmT3vrkSVM4crSY46e3kj68oIeo20cNe/jwaDGdnVeYOXU5ifHD6LjU3FWMFcrXG2xCMt3RhV4tczNyzkMpj90b4aj9brfbmfqZ6cTYb/a4Gt+y6xekJk/RrSmSOOw0paXFFBYW0XI2hcy0gi6XjhbN8uqaV3j//eqQuVee6H59Wjlct5WkYQ1sLN1g4sgcaKkfQFwifXQ677/3Hm3trcTHJTB9eg4FM3PZurWSnTt2dL0eF59ARsqdIXu9ZhLW6Y4u9AY8zWhlFi5l5WCe9ruZWCwWVv30x/22W/PXp+oK1NntnZrT5LLSZzJl3H1hkSbnS7u6QLWo06q3np5cSE1NDe/ufYfOThvNLRco3VTMD37wAzZtKqG55ULX6y+//GLIXm+oENIrdr2YKbPqy2rYX8nhaMFbd6Q//GMx99/1IhaLVccxbfxp/Tfp7LR1HX+QuIkz5w6xsODHmt1qm3f/hBdfei4o1Zla0dtd6shH5R5b1IXDuEL1egNBQFfsQogFQohDQogjQojvG3FMXzAz51zvajgURMPChe7dkY581HclNmBAnF8CV67jHz9VpTtNbtSNedx99xdCWoPE2/1z4QpE1tS9yYYNb5lu5MyqHwjV6w0l/F6xCyGsQC0wDzgJvAN8WUr5oad9zFqx+6ql8h1LBbbO/rv96CGYPVbDmZ6+2PyuQNqWXb/QrSnizqeaEJ/EooKndGuWv/H2f/JPRc+HvAaJp/t3+Uo79Y3vUtdYQVyChQ0b3grIuAvnFdFyTr/v3xUbcUf3IOz2bdtoa2/FahnANUkpTMi4k1HDp3DFdiko1xsItK7YjciKmQockVIec574FeCzgEfDbhah0Mqst+Swp3Z8Lgld2Qx3FRYFJdgaarhvt9ZKzMBBtHWc9xjs7I3Lp/q71c/3eL2jo82nNDlbZ0cPn3DGqHyO1lcwfXpuSGmQeLp/8XEJTM/J4cXV7lvUmUXVju0szFuta5/U5NvYUPma2/d6B2EXFTzVIwj7Xs3r7Nj3W2JiBjAjNy/g1xtKGHHFKUB9t79POl/rgRDiW0KIPUKIPWfPnjXgtH0xo2JVL0YHcIPd3CPQuOQISkuLuwJmrW3NDBmawNF6bZoiR+sriE+09snTd6UF6qGvTGxoa5C4u3/NLRcoLS32WIRnFkbm12sJwn52ztNMveWrxMbF89xzawN+vaGEEVftbgnVx78jpXxBSnm7lPL2YcOGGXDavhhRseovRoqGKT+9A6N8qjk5M6hv3Kvr3PWN73LD9WPdvheNGiR68HUi7Z1vHukibmZghGE/CaR2+3skcNqA4+omFFqZGRXA7d3cI5fhnKCZ59nP99nJ76nhdHMTV+rOMXXybdTU1Bh2DaFIVlYWO3Zso/FCBZt3/8Stpsjm3T/hzMVKj+4RX9ICa46/zdjRc9y+H40aJHrwdSLt3fwj0kXczMAIw/4OkCmEGC2EiAHuA9404Li6CYWccyNEw3r76c/Qzkp28TpHmcRQVjONF5jJaqYxkxQS2yS3T7w1Koz74cM1vPjScyQOO82Gykf50/pvsqHyURKHnebFl56jtvagR593d30SLRz5yJH7PmLYRI/bRJsGiR6Myq+PdBE3M/DbokkpbcBSYCNwEPirlPKAv8f1lWC3MjNCNKy7n75RtvEUe5nPKB7nM27b8j3BVO7pHM3026dGhVvGVx+yHpfO4bqt7Dv4OrOm/juinyyraNMgcaGlxaDeidRTbCQcRc6CjSFaMVLKDUDI1Oq601Lp3unn+UeeNS1aboRomMtPL4Ff8T6f56Z+i66EEBSQgmwVKsPGCy6XTlHRIk40lPdJC/yoYQ+1J7bQabcxP/cxBif2LyMRbRok0Dc7ZWHe6h7ZKUseWNaVDtpb6EtLp6Pen91wFDkLNiEpAmYErpVdoCsGjRAN276zilVM4QBNxGAlD+8aNQAFDGd704Eu8SuFezylBVotA0hKGM7k8V9ixLCJ/a7UXQSzIXQw6NFicGJPQ+0pHbS/ibR7vrmn2Eg4iJyFGmpZZzBGBHBdfvoyTjGLlJBoyxdpuHPpvPnW34mNG8CIYTdrMurRpkHia3ZKRkaGX7ERo4Kw0YQy7AZjRADX5ac/xKch05YvGjDKJ6wHLb7qUMGf7BR/YiOhLHIWqkSUCFgo4Y+ErkvM7PfU8AIzgy6REE30cDVo9An7GoDXImcbStIFZkgEaCGaRb96o1VSQBl2E7Hb7V0B3KpeAdwHH/mexwBuSUkJS+9dzOnmJlYzjSShXSKhWV5mZVw1F9pajLyUqCIQmiuBnECMIjFxMAvzVuvW2tlQ+SjNLRf8Onc43i8ziAg99lBDb3m/6/GzeMsmLrS1aOqxClf99DcSF3SJhGjE33x5b4RrJWUws1OMKFALFsFwt0VsVozRGNlH1RsuP/3Uybexqa2ePBk6bfmiBTOzqrp81RN9a/UXLIKdnRJqImda0JMaauRkpFwxGgiWDG9NTQ23TbyVL3WOpkD00VXrQzDb8im0Eyxftb/40oKv9kQZew/+hYKCfB56aHlUNZcxw32kXDEG0bu8v7/H5jw5nIXNydxVWGTI41V2djbv7n+PfyScZqs8HdJt+RTaCddKSl+yUw4d38y0mx+g5WwKSx5YRmZmdsRXR0Pw3W3KAnjBjD6qesjOzmbX3j3sSLsUFIkEhfGEayWlT1o7spO0lM909TZNviaf6dNzI964B1u4TBl2Lxgpw+sr4dikWuEZo+RsA42/Wjuu1WlW6iImTrgVi7CEbM6+vwRbuEz52L0wOD6BVe1TVMqhwjB88VW7a/UXLLylgx46vplOu41ZU//drdaOlJJ1W1cyedw9DL1uTEjm7PuLWamhgWyNF9EYIcOrUHTnoYeWs+SBZX63+gsWvbNT3nj7j9hsVxg4MJYbrh/rVWtHCEH2TXOprStjZPKkkG436CvBdrcpw+6FUOijqogsuvuqtVRSGiFdYDTd00ETEwfz+bnP6FqdpibfxrsHXun62+WmAbjzzrvCvmo02Kmh4XvnAkQo9FFVRBZGtfoLFXxdnV6xdfR5PVK6HwVbuCw0PykhRCj0UQ1Xoq0Rtx7CuZKyN0Y1CYfI6X4UbOEyZdi9EAp9VI0mEAZXNeL2jtnSBYHC6Cbhvubsh5JSZjCUQrujsmI04Ko8XdicTF4/laeVlkbWG1h5aga9pREmM9QhjYCNas5RmdTUQxrBbrdTWlrK2qfXsH1nFa0d7STExjFjWg7LH17htpIwWJW6iuBQUlLCkgeWMWfqSs3B4HXlK5ky/l5Sbrilz/t2u40/rf8mnZ3aW0yGolJmMCtPlWHXiD8yvKGCXoP7h7/8mWXf+o7mSQAcq6bsmzLIrQqfwEsAAAw1SURBVI8lT3rv/KRkEMIX16S/Zs2zlG8t57YJXyYrfZbX/Q7XbeXDoxu5e9ZP3WbO6FWEDGXlR6OVQpVhNwFfZXhDAb0G9y1OUMxHfFlk6lp1l5SUsOzexTzaPEHz6u3JpAM8/+rLqp1fGNF7hXxtUgplu9cyKfsLZKbN9Ph5OfJROfsOvt5vP1k9OfvhoNXushtr1jzLjqqqHsJlK1Y8qMtuqDx2EwhWH1UjuCqNkA5e7K1dSnbSyD8xxtFv1cP2Ln0c2UxXE21/KnXD8b5GI576ni7IfYyyXWupPbGFsaPn+NQkXG/OfjgoZQbDbqgVe5Tg6sqUL0Z43fYDeZ6/cYwfcrvuVfeXPv8FVakbwXhbIUtp5/TZ/Rw6vpmPzx/iypV2hMXKdYNHMmncPV6bhOtdUYerUqavqBW7ogfbd1axCm2Kgv400VaVupGNtxWyEBZSbrilKygqpWRj1Uo6Ll+gveMTPD3+9faBa3VNVO3YzsK81bquITX5NjZUvqZrn3AjNB3CCsPRY3D9aaLtasSth3Ct1A2l9LpA4Yu41ZiRc5k4caIpOfvBLt0PVdSKPUrQI43Qgc3nVfe8gtlUl50kH+8uHxfhWKkbrM44wcafFfKFi58Y3v1Ia+m+lHZOf7yfQyc2c+ZcDXZ7J4mJg8nJmRGRDUCUYY8SZkzL0WxwY50pjXr1cYQdNpa9zQ3ERXQ7P0/BQ4DYQUkRI2rlSmd85pm1VFVtp72tFYSF7dX/Q/bouYy4oX9/uQvXCtmMIKKrOKo/H/vFlgbKdq3Fah3I2NFzyZn0jYifhCNnilL0ix5phLFcq1sfZy9nyeZa/h8FCKCS05r2C4dK3e4EuzNOoKitrSUzM5slDyyj5WwKC/NW85VFv+WewrWMGn4b1Qdf5c0tj3GxpcHrsczUkvdWun+xpYGN255kfEYRCwueIDOtgNhBSVgs1q5JOBIbgCjDHiXokUaYRQqbqNelc1HGKQpJZaCw8u/cyhscp1yeirh2fsHujBMIXE8kydfkM2fqSrfGcGHBE4zPKGLjtie9Gncjxa1601/pvpR2ynatZdK4e8hMi9xJ2B3h8W1S+I3FYmFdaQnrkhqpEA39GtwmcYnz4hKVWvVxaMCGZDxDAEgW8TzCFEqp58e8E1Ht/ILdGcds9DyRZKYVMGncFynb/UukdG8MjRa36k1/SpmnP96P1RpDxqjInYQ9oQx7FJGVlcW23TvZntrRb//UHaMu8/fidazXMAlUyNO8wTGWcjOWbkYgWcTzE+7gHsZQzVm+RxXfFuVh384vXBtRa0XvE0nGqAKslgGcPrvf7fv9iVsZlVXkSSmz5vjbjB09J2In4f5QBUpRiFZphP70cfZyljJOYUOylJtJFv2nnFVwmsaZIyneEt6rIavFylcW/RaLxap5H19ErYKFrwU/9Y3VzL7jP7pe86bNYoZoV9/S/TbuKVxreHu6YKIKlBQe0Zqd4Gqi7ZoEVjonAWGHbK7lHsYwniE9VuqemCyHsnLnDqMuIWgEuzOO2fiazrhn/yvY7bY+4laejLoZWUW9P9dWizVqc9yVYVf0i7tJwGqxsFzeglVDqpuLSKku1ZJe1xszg4dG43s3pHb+tP6bXvPSe/vwPWFEq7xIn4T7wy8fuxDiS0KIA0IIuxDC6+OBIjKIpurS3gS7M47Z+NoNKTEhic5OG80tFygtLWb+/PluDXEgs4qC3Z4umPgbPN0PfAHQ1iZEERFEcx9YvZ1xak+U0dxyjs7OzrBIozPbGAYyqyjSJ+H+8MuwSykPSikPGTUYRXgQzX1g9TSirj1RxnuH/kZG6jy+uWQ5mZnZIV8AY7YxDGRWUbDb0wWTgPnYhRDfAr4FMGrUqECdVmECXcVOrY2amnaEW3WpN1zpdUVFizjRUO62M86h45vptNtYkPsDBicOR8q7wkJioLsx1NK4Qq8xDKRol2sSnj49F0Bzd6VwKZbrD6+GXQjxNpDs5q0fSCn/ofVEUsoXgBfAke6oeYSKkMNV7JQ7dRqyGW19YEt3RsQXxoWrEfXGjRv5p3vv5539f6Kz8zIDB8Ryw/VjmTz+Sz20x40IBgYCs41hoAOaWibh/jJ4whWvhl1KOTcQA1GEF65ip0WFC6hsOuC+D2ySow/sttLwqi7VisViQQjB4MTr+dxsbY2cg9HBRy9mGsNgZBW5JmGjlSVDGUMKlIQQW4HvSSk1VR2pAqXIIZz7wBpBJHfwMbJXp4uSkhKWPLCMOVO1TYRSSt7e9QS/+/3zITsRBpKANLMWQnwe+BUwDPgU2Cel9Hr3lWFXRAqJiYNZmLc6oqobzSQcmk+HMgGpPJVSvgG84c8xFIpwRnXw0Uc0BzQDiao8VSj8IJqrG30lWgOagUQZdoXCDyJdYsAsojGgGUiUYVco/OChh5az5IFlZIxy71Lojaug53ernw/A6EIbM1rlKRyo6VCh8INorm5UhC5qxa5Q+IEKBipCEWXYFQo/UcFARaihDLtCYQAqGKgIJVRrPIVCoQgTtBYoqeWDwbga9BbNnsfg+ASsFguD4xMomj1PV4NehUKh8BXlijGQ2tpaFs1bgP0TR/PnVUwhngG0tduoLjvJ0j2LsQxJYF1pifKzKhQK01CG3SBqa2vJnTqNRc3J5NrTe2RGJBFDPiPIax7OttZGcqdOY9vuyFQ8VCgUwUe5YgzAbrezqHABi5qTyZPDPRaqCCHIk8NZ2JzMXYVFyi2jUChMQRl2AygtLUU2tZJrd9ePpC959mQ6m1p8atCrUCgU3lCG3QDWPr2G3OYhuhr05rUM4dmnfm7yyBQKRTSiDLsBbN9ZxWSG6tpnshxK1c4dJo1IoVBEM8qwG0BrRzvxOuPQcQyg9VK7SSNSKBTRjDLsBpAQG0cbNl37tGMjYVCcSSNSKBTRjDLsBjBjWg7VnNO1T7U4R8606SaNSKFQRDPKsBvA8odXUJnUhFZ5BiklFQlNPPjI90wemUKhiEaUYTeAwsJCLEMS2GZp1LR9paWRAdcnKk1uhUJhCsqwG4DFYmFdaQnrkhqpEA0eV+5SSipEA+uTGnmrtFgp/SkUClNQkgIGkZWVxbbdO1lUuIDKpgPkNQ9hMkOJYwDt2KgW56hMasI6JJFtpUpOQKFQmIdaMhpIVlYWNceO8PyrL9M4ayQr46r5jqWClXHVNM4cyfOvvszBY4eVUVcoFKai9NgVCoUiTFB67AqFQhGlBGXFLoQ4C9QF/MQ9GQo6k88jF3UvrqLuxVXUvbhKqNyLNCnlMG8bBcWwhwJCiD1aHmmiAXUvrqLuxVXUvbhKuN0L5YpRKBSKCEMZdoVCoYgwotmwvxDsAYQQ6l5cRd2Lq6h7cZWwuhdR62NXKBSKSCWaV+wKhUIRkSjDrlAoFBFGVBt2IcSXhBAHhBB2IUTYpDIZiRBigRDikBDiiBDi+8EeT7AQQrwohPhYCLE/2GMJNkKIVCFEmRDioPP7sTzYYwoWQohYIcRuIcR7znvx42CPSQtRbdiB/cAXgIpgDyQYCCGswPNAETAe+LIQYnxwRxU0XgIWBHsQIYINWCGlHAdMA74bxZ+LS8BsKeWtwCRggRBiWpDH5JWoNuxSyoNSykPBHkcQmQockVIek1JeBl4BPhvkMQUFKWUF0BTscYQCUsoGKeVe5+/NwEEgJbijCg7SQYvzz4HOn5DPOIlqw64gBajv9vdJovQLrHCPECIdmAzsCu5IgocQwiqE2Ad8DGySUob8vYh4PXYhxNtAspu3fiCl/EegxxNiCDevhfxqRBEYhBCJwOvAg1LKi8EeT7CQUnYCk4QQ1wJvCCEmSilDOhYT8YZdSjk32GMIYU4Cqd3+HgmcDtJYFCGEEGIgDqP+Rynl34I9nlBASvmpEGIrjlhMSBt25YqJbt4BMoUQo4UQMcB9wJtBHpMiyAghBPBb4KCU8plgjyeYCCGGOVfqCCHigLlATXBH5Z2oNuxCiM8LIU4C04H1QoiNwR5TIJFS2oClwEYcAbK/SikPBHdUwUEI8WdgBzBWCHFSCLEk2GMKIjOArwGzhRD7nD93BntQQWI4UCaEeB/HQmiTlHJdkMfkFSUpoFAoFBFGVK/YFQqFIhJRhl2hUCgiDGXYFQqFIsJQhl2hUCgiDGXYFQqFIsJQhl2hUCgiDGXYFQqFIsL4/wFhkS8dmXS91AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 导入matplotlib绘图库\n",
    "import matplotlib.pyplot as plt\n",
    "# 导入生成分类数据函数\n",
    "from sklearn.datasets.samples_generator import make_classification\n",
    "# 生成100*2的模拟二分类数据集\n",
    "X, labels = make_classification(\n",
    "    n_samples=100,\n",
    "    n_features=2,\n",
    "    n_redundant=0,\n",
    "    n_informative=2,\n",
    "    random_state=1,\n",
    "    n_clusters_per_class=2)\n",
    "# 设置随机数种子\n",
    "rng = np.random.RandomState(2)\n",
    "# 对生成的特征数据添加一组均匀分布噪声\n",
    "X += 2 * rng.uniform(size=X.shape)\n",
    "# 标签类别数\n",
    "unique_lables = set(labels)\n",
    "# 根据标签类别数设置颜色\n",
    "colors = plt.cm.Spectral(np.linspace(0,1,len(unique_lables)))\n",
    "# 绘制模拟数据的散点图\n",
    "for k,col in zip(unique_lables, colors):\n",
    "    x_k=X[labels==k]\n",
    "    plt.plot(x_k[:,0],x_k[:,1],'o',markerfacecolor=col,markeredgecolor=\"k\",\n",
    "             markersize=14)\n",
    "plt.title('Simulated binary data set')\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100, 2) (100,)\n"
     ]
    }
   ],
   "source": [
    "print(X.shape, labels.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100, 3)\n"
     ]
    }
   ],
   "source": [
    "labels = labels.reshape((-1, 1))\n",
    "data = np.concatenate((X, labels), axis=1)\n",
    "print(data.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X_train= (90, 2)\n",
      "X_test= (10, 2)\n",
      "y_train= (90, 1)\n",
      "y_test= (10, 1)\n"
     ]
    }
   ],
   "source": [
    "# 训练集与测试集的简单划分\n",
    "offset = int(X.shape[0] * 0.9)\n",
    "X_train, y_train = X[:offset], labels[:offset]\n",
    "X_test, y_test = X[offset:], labels[offset:]\n",
    "y_train = y_train.reshape((-1,1))\n",
    "y_test = y_test.reshape((-1,1))\n",
    "\n",
    "print('X_train=', X_train.shape)\n",
    "print('X_test=', X_test.shape)\n",
    "print('y_train=', y_train.shape)\n",
    "print('y_test=', y_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch 0 cost 0.693147\n",
      "epoch 100 cost 0.554066\n",
      "epoch 200 cost 0.480953\n",
      "epoch 300 cost 0.434738\n",
      "epoch 400 cost 0.402395\n",
      "epoch 500 cost 0.378275\n",
      "epoch 600 cost 0.359468\n",
      "epoch 700 cost 0.344313\n",
      "epoch 800 cost 0.331783\n",
      "epoch 900 cost 0.321216\n"
     ]
    }
   ],
   "source": [
    "cost_list, params, grads = logistic_train(X_train, y_train, 0.01, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'W': array([[ 1.55740577],\n",
       "        [-0.46456883]]), 'b': -0.5944518853151362}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.]\n",
      " [1.]\n",
      " [1.]\n",
      " [0.]\n",
      " [1.]\n",
      " [1.]\n",
      " [0.]\n",
      " [0.]\n",
      " [1.]\n",
      " [0.]]\n"
     ]
    }
   ],
   "source": [
    "y_pred = predict(X_test, params)\n",
    "print(y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.27936388],\n",
       "       [0.92347241],\n",
       "       [0.90814155],\n",
       "       [0.24192851],\n",
       "       [0.94789076],\n",
       "       [0.98545445],\n",
       "       [0.29242059],\n",
       "       [0.06128463],\n",
       "       [0.95821867],\n",
       "       [0.46214622]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigmoid(np.dot(X_test, params['W']) + params['b'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       1.00      1.00      1.00         5\n",
      "           1       1.00      1.00      1.00         5\n",
      "\n",
      "   micro avg       1.00      1.00      1.00        10\n",
      "   macro avg       1.00      1.00      1.00        10\n",
      "weighted avg       1.00      1.00      1.00        10\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import accuracy_score, classification_report\n",
    "# print(accuracy_score(y_test, y_pred))\n",
    "print(classification_report(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n"
     ]
    }
   ],
   "source": [
    "def accuracy(y_test, y_pred):\n",
    "    correct_count = 0\n",
    "    for i in range(len(y_test)):\n",
    "        for j in range(len(y_pred)):\n",
    "            if y_test[i] == y_pred[j] and i == j:\n",
    "                correct_count +=1\n",
    "            \n",
    "    accuracy_score = correct_count / len(y_test)\n",
    "    return accuracy_score\n",
    "\n",
    "accuracy_score_test = accuracy(y_test, y_prediction)\n",
    "print(accuracy_score_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8888888888888888\n"
     ]
    }
   ],
   "source": [
    "accuracy_score_train = accuracy(y_train, y_train_pred)\n",
    "print(accuracy_score_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.31055924594920103"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train[1][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.17221697, -0.73380144],\n",
       "       [ 2.54116921,  0.31055925],\n",
       "       [-0.00718884, -0.70554359],\n",
       "       [-0.31285288, -0.17275221],\n",
       "       [-0.67290531,  0.79310561]])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd4VGXax/HvQw8BRAWCIBCKIoJJ1FDEihVbYgdX3bUtu+6ugmJ3QSXrawcsq7u4KrurSxVQULEgKootSDIISO+Q0FsKac/7xyRAwiSZyZQz5fe5rrniJCcz9xniuc/T7sdYaxEREanndAAiIhIelBBERARQQhARkXJKCCIiAighiIhIOSUEEREBlBBERKScEoKIiABKCCIiUq6B0wH4olWrVjYxMdHpMEREIsqCBQu2W2tb13ZcRCWExMREMjMznQ5DRCSiGGPWeXOcuoxERARQQhARkXJKCCIiAighiIhIOSUEEREBlBBERKScEoKIiABKCCIifpn9yxam/bzR6TACQglBRKQO9h8o4YEp2fzxnZ+Z9NMGomF/+ohaqSwiEg5+Xr+LeydlsWFnPn8Z0I2hF56AMcbpsPymhCAi4qWS0jJenbuSV75YSdsWTZg45Az6dD7G6bACxtGEYIy5F7gTsMAi4DZrbaGTMYmIeLJuRx7DJmWxcP1urj61PU+m96RFk4ZOhxVQjiUEY0x74B7gZGttgTFmMjAYGO9UTCIiVVlrmbpgI098sJj69Qwv33gqacntnA4rKJzuMmoAxBljioGmwGaH4xEROWhXXhGPTl/Ex7/k0LfzMYwelEL7lnFOhxU0jiUEa+0mY8wLwHqgAPjUWvtp1eOMMUOAIQAdO3YMbZAiErO+WbGd4VOy2JlXxMOXnsTvz+5C/XqRP3BcE8emnRpjjgbSgc5AOyDeGHNz1eOsteOstanW2tTWrWvd30FExC+FxaVkzFrCzW/+QLPGDZj+pzP547ldoz4ZgLNdRhcCa6y12wCMMdOA/sA7DsYkIjFsWc4+hk5cyK85+/jtGZ145NIexDWq73RYIeNkQlgP9DPGNMXdZXQBoO3QRCTkysos4+ev5ZnZv9KiSQPevrU3A05q43RYIefkGMIPxpipwM9ACbAQGOdUPCISm7buLWT4lGzmrdjOBSe14dnrkmjVrLHTYTnC0VlG1trHgcedjEFEYtfsX3J4ZJqLguJS/nZVL27q2zEqVhzXldPTTkVEQi7vQAmjZi5hUuYGerVvwdhBp9KtTTOnw3KcEoKIxJSF5XWI1u3M50/ndWXYhSfSqIHqfIISgojEiJLSMl77chUvzVlBQvPGTPh9P/p1OdbpsMKKEoKIRL0NO/MZNimLBet2kZ7SjlHpvTgqLrrqEAWCEoKIRC1rLe/9vIknPliMAV4anEJ6SnunwwpbSggiEpV25xfx2Ixf+NC1hT6JxzB6UDLHH93U6bDCmhKCiESd+Su3c9/kbLbvP8ADl3SPmdIT/lJCEJGocaCklBc/Xc4b81bTuVU80397Jqccf5TTYUUMJQQRiQrLc/cxdGIWS7fs5aa+HXns8h40baRLnC/0aYlIRLPW8u/5a3n6419p1rgB//ptKheenOB0WBFJCUFEItbWfYU8MMXFV8u3MaB7a567LpnWzWOzDlEgKCGISET6dHEOD09b5C5Dkd6TW/p1iuk6RIGghCAiESW/qISMWUuY8OMGerZrwUuDU+jWprnTYUUFJQQRiRjZG3YzbFIWa3fk8YdzuzD8ou6qQxRASggiEvZKyyyvf7mSsZ+voHXzxvzvzn6c0VV1iAJNCUFEwtqGnfncOymLzHW7uCLpOJ666hSOaqo6RMGghCAiYclay/SFmxj5vrsO0ZhByVyV0l4Dx0HkaEIwxrQE/gX0Aixwu7X2OydjEhHn7ckv5rEZi5jl2kLvxKMZfUMKHY5RHaJgc7qF8BIw21p7nTGmEaB/cZEYN3/VdoZPzmbbPtUhCjXHEoIxpgVwDnArgLW2CChyKh4RcdaBklJGf7accV+vJvHYeN67qz/JHVo6HVZMcbKF0AXYBrxtjEkGFgBDrbV5DsYkIg5YuXUf90zIYsmWvdzYpwMjrjhZdYgc4OQE3gbAacDr1tpTgTzg4aoHGWOGGGMyjTGZ27ZtC3WMIhJE1lr+891aLn/5G3L2FjLultN5+pokJQOHOPmpbwQ2Wmt/KH8+FQ8JwVo7DhgHkJqaakMXnogE07Z9B3hwajZzl23j3BNb8/z1SbRp3sTpsGKaYwnBWptjjNlgjOlurV0GXAAscSoeEQmdz5fk8tB7LvYfKOHJtJ789gzVIQoHTrfL7gbeLZ9htBq4zeF4RCSI8otK+NuHS/nfD+vpcVwLJgxO4cQE1SEKF44mBGttFpDqZAwiEhqLNu5h6MSFrNmRx5BzujD84hNp3KC+02HJYZxuIYhIlCsts/zjq1WM+Ww5rZo15t07+tK/WyunwxIPlBBEJGg27srnvknZ/Lh2J5efchxPXd2Llk0bOR2WVEMJQUSC4v2sTfx1+i9Y4MXrk7nmNNUhCndKCCISUHsKihkx4xc+yN7M6Z2OZswNKXQ8VlVpIoESgogEzA+rd3Df5Gxy9hZy30Un8qfzutKgvjawiRRKCCLit6KSMsZ8vpx/fLWKTsc05b27+pOiOkQRRwlBRPyyatt+hk3MYtGmPQzu7a5DFN9Yl5ZIpH81EakTay3v/rCev324hLiG9fnHzaczsFdbp8MSPyghiIjPtu8/wMPvufh86VbOObE1L1yXRJsWqkMU6ZQQRMQnc3/dygNTs9lbWMLjV57M785IpJ42sIkKSggi4pWColL+76Ol/Pf7dZzUtjnv3tmP7m1VhyiaKCGISK1+2bSHYZOyWLl1P3ee1ZkHBnZXHaIopIQgItUqLbO8MW81L366jGPjG/POHX056wTVIYpWSggi4tHm3QXcNzmL71fv5LJT2vJ/V5+iOkRRTglBRI7wQfZm/jp9EaVlluevS+K6049XHaIYoIQgIgftLSzm8fcXM33hJk7r2JIxg1LodGy802FJiCghiAgAP63dybCJWeTsLWTYhSfwlwHdVIcoxighiMS44tIyxn6+nNe/XEWHY5oy5Y9ncFrHo50OSxzgeEIwxtQHMoFN1tornI5HJJas3rafeydlkb1xDzekHs/IK3vSTHWIYlY4/MsPBZYCLZwORCRWWGuZ8OMGMmYtoXHDevzj5tMY2Os4p8MShzmaEIwxxwOXA08B9zkZi0is2LH/AA9PW8RnS3I5+4RWvHB9MgmqQyQ430IYCzwIaP27SAh8uWwrD0x1sSe/mL9e3oPbz+ysOkRykGMJwRhzBbDVWrvAGHNeDccNAYYAdOzYMUTRiUSXwuJSnv5oKf/+bh3dE5rzn9v70OM49dJKZU62EM4E0owxlwFNgBbGmHestTcffpC1dhwwDiA1NdWGPkyRyLZ48x6GTcxixdb93HZmIg8NPIkmDVWHSI7kWEKw1j4CPAJQ3kK4v2oyEJG6Kyuz/Oub1Tz/yTKObtqI/9zeh3NObO10WBLGnB5DEJEg2Ly7gOGTs/lu9Q4u6ZnA09ckcUy86hBJzcIiIVhrvwS+dDgMkagwy7WZR6ctoqTM8uy1p3BDagfVIRKvhEVCEBH/7SuvQzRt4SaSO7Rk7KAUOrdSHSLxnhKCSBTIXLuTeydnsWlXAfec3427LziBhqpDJD7SX4xEHFeui/SJ6XQe25n0iem4cl1Oh+SY4tIyRn+6jBv++R0AU/54Bvdd3F3JINRcLkhPh86d3V9dkfk3qRaCRBRXrov+b/Ynvzgfi2XdnnXMWT2H+XfMJykhyenwQmrN9jyGTcoie8Nurj3teJ5IO5nmTRo6HVbscbmgf3/IzwdrYd06mDMH5s+HpMj6m9RthNSZE3fqI+aOOJgMACyW/OJ8Rs4dGfT3DhfWWib+uJ7LX57H2u15/P03p/HiDclKBk4ZMeJQMgD31/x8GBl5f5NqIUidhPJO3ZXrYsTcEbhyXOTm5R5MBhUsluzc7IC+Z7jamVfEw++5+HRJLv27HsuLNyRz3FFxTocV21yuQ8mggrWQHXl/k2ohSJ2E6k69IvHMXDaTtXvWUlBScMQxBkNyQnJA3zccfb18G5eM/Zovl23jsct68M4dfT0mA42xhFhSElSd1msMJEfe36RaCFInrhxXSO7UqyaeqgyGpg2bMmrAqIC+bzgpLC7lmY9/Zfz8tZzQphn/vq0PJ7fzXIdIYywOyMhwjxlUdBsZA02bwqjI+5tUC0HqJKltEobKd0XBuFP3lHgAmjZoSmLLRNK6p1W62AXr7tipu+6lW/aS9uo3jJ+/llv7JzLz7rOqTQagMRZHJCW5B5DT0iAx0f01AgeUQS0EqaOMARnMWT3n4MUnWHfqSW2TWLdnXaWkYDBc1PUiZgyeUenYYN0dO3HXXVZmeevbNTw3exkt4hoy/rbenNe9Te2xhqjlJlUkJcGMGbUfF+bUQpA6SUpIYv4d80nrnubxTj1QMgZk0LRh04OtkZoST7DujkN9152zp5Bb3vqBv324lHO7t+aTYWd7lQwgdC03iU5qIUidJSUkHXGXHoz3mH/HfEbOHUl2bjbJCcmMGjDKY+IJ1t1xKO+6P1q0hUemLaKopIynrzmFwb19q0MUqpabRCclBAl73iae6rqX/L07DtbrHm7/gRKe+GAxUxdsJPn4oxgzKIUurZv5HqsPCVSkKmOrzp8NY6mpqTYzM9PpMKLO4fP8k9omkTEgIyIvIFX7+ivujgM9hhCo162wYN0u7p2UxcZd+fzpvG4MvVB1iCSwjDELrLWptR6nhBDbgn2xCzVXruvg3XHnlp0ps2Ws273O70R3+OsG6q67pLSMV75YyatzV3LcUU0YMyiF3onH+PWaIp4oIYhXd/7pE9OZuWzmEd0had3Tgj4+EEzhnujWbs/j3slZLFy/m2tObc8T6T1pEe6lJ1wud5kGl8s9qyYjIyKnVsYibxOC2qVRquoK35nLZtL/zf5HzJ+P1mmK4Tof31rL5J82cNnL81i1dT+v3Hgqo3vUo8Wg646olBlWK44rCrjNnAlr17q/9u8fsVU9xTMNKkepEXNHkFecd/C5xZJXnMfIuSMr3fmHYsDUCeGY6HblFfHItEXMXpxDvy7HMPqGFNqtX+GxUqbr4/H0n3dr+Kw4rqmAWxTMvxc3x1oIxpgOxpi5xpilxpjFxpihTsUSjTI3ee5a+2nzT5We+zLPP5L4PR8/wPXt563YxsCXvmbOr7k8culJ/O/OfrRrGVfthXbEtLvDq4UTRQXcpHpOdhmVAMOttT2AfsCfjTEnOxhPVCm1pR6/X2bLKj0P1QKz2gS6e8SvRBfA7pHC4lIyZi3hljd/pHmThky/qA1/GH0f9bp2cSeazEyPF1pX/W3h1cLxp4BblGweEwsc6zKy1m4BtpT/9z5jzFKgPbDEqZiiSX1T3+P365kj7wFCscCsJsEoDVFpPv6Gn0jeXMaor+qR9MOI2gdDA9Q98mvOXoZNzOLXnH387oxOPNyhhLhzzqrcPVSvnvvCenhSMIak0tasIzd8uvIyMuCzz6DgsGqzTZrUXsAtijaPiQVhMahsjEkETgV+cDaS6JHaPtVjl0nvdr0diqh6wRoATkpIYsbJo1jz5B5mvJRL0sLN3t3t+9k9UlZmefObNaS9+i3b9xfx9q29eTK9F3FPPn5koiktPZQU4GClzIxrXomOrrwo2jwmFjieEIwxzYD3gGHW2r0efj7EGJNpjMnctm1b6AOMUJE0NhDUAeC6XJD86B7J3VvI797+kYxZSzjnhFbMHnY2A04qr0PkKdEAJCQcqpR53nmQmkrSbx9gviuVtOPOc7Qr76ARI6CwsPL3Cgtrv7Br7CGiODrLyBjTEHcyeNdaO83TMdbaccA4cK9DCGF4ES2SShgEdaZTXS5IdaxvP/uXLTw8bRGFxaU8dXUvftOnY+U6RElJ7i6TKt1D9O7t7oqq0r2StG4dM2Y1rb57JZTrAup6Ya/unCNw85iYYK115AEY4D/AWG9/5/TTT7cSfbJzsm38U/HWPGEsT2DNE8bGPxVvs3Oy/X/xtDRrjbHWfUlyP4yxNj29lqCy3cckJrq/Zlcfy/7CYvvAlCzb6aFZ9oqX59mVW/dV/5rx8YfiMcb9vOK1fYm1ttcKNH8+x1DGKR4Bmdab67I3BwXjAZwFWMAFZJU/Lqvpd5QQold2TrZNn5BuE8cm2vQJ6YFJBtYG/YL087qd9pznvrCJD8+yz3681B4oLq09nuoSTWJi5QtuxSMx8cjXqesFuq78+Rx9SK5ex5KW5n69tLSaX8+XY6NY2CeEujyUEGJHdk62TZuQZhPHJNq0CWn+JYhAX5CstcUlpXbsZ8ttl0c+tP2fnmN/WL3D79f06SLvS/LwlaeLaHa2teeea21cnPsxYIAzF1dfEpNaJwcpIUjE8rcLKaDJxIN12/Ps1X//xnZ6aJYdOuFnu6egKDAv7MsFLFgtBE8xVCSBcLiw+nLeoW5FhTFvE4Ljs4wkvIRD/Rx/pqHWVsPJn/Oz1jJ1wUYufelrVmzdz0uDUxg7+NTAFaXzZW/ejAz3QHeV6ao+bezuacGYp1lZBQXux+Hfc2rqqC+D25rh5DPVMpKDnNg72GMcfkxDrSmZjBowqs7ntzu/iEenL+KjRTn0zdvM6C9ep/3CDoGf2ePt3rwVyWPkSPcFLjnZnQy8jaW6BWNHHeV5amxVTl1YfZm1pBlOPquxhWCMaWGM6erh++E3dzGGBOsu3skKoYefU0FJwRGL6gBy9+fWer41JZO6nt+3K7czcOw8Plucw4Pz/8f/XruL9ksWOl/xsyJ5rFnj/upLYqpufUZZ2ZFrMDxx6sLqS8soEK2oGFNtQjDG3AD8CrxXXnzu8CWu44MdmHjmbVnrOr12He7MA5Gcqp5Tbl7uwT0MDldQUlDr+dZU1M7X8ztQUspTHy7hpn/9QNPG9Zm+6j3+9M0E6peV14lysuvEX9V1p9Sr575o1sTJC6sv3Wq+HCtAzS2ER4HTrbUpwG3Af40x15T/zPtdvyWggnkX72uF0EAlp6rnVCGhWQJxDeIqfa+2861phbYv57c8dx/pr37LG/PWcEu/Tnx499n0+nFu9PRJV7cau3dv90WzbVvPvxcX5/nCGsoCdr60jPxpRcWgmhJCA+suQIe19kdgAPCYMeYewItORgmGYJZ58LXcRV2TU9VWReamzCPOCaBJgyYkxCcc8f2azrem6q3enJ+1lvHfruGKV75h+/4DvHVrKhlX9SKuUX3/Kn6Gm5q6U5KS4JNPID6+8s/j4+H774+8sHpTHVYVTyNCTYPKe40xXa21q8BdndQYcx4wA+gZiuDkSMEs8+BruYu6djFVHdj1VIG14pwqjvHlfKur3lrb+W3dW8j9U118vXwb55/UhmevTaJ188aHXqCOJS3CUm2D0r4MWtdWHVYVTyNGtXsqG2POBrZYa1dW+X5D4BFrbcj/L9Ceyu4Lar9/9aOg5FAZ4rgGcXx/5/chr1NUl/2Yq/udeqYeZbbsiP2PgZDsjfzp4hwenraI/KISHrv8ZG7uW6UOUQWXq+4ze6JV587ulkFViYnurpr0dHeroepsn7Q07bYWIoHYU/nfwLXGmIOtCGNMAvA2cKX/IUq48XWA+JakWyrd3XtTUbW6VkVCswSP3TzB3sAnv6iER6a5GPLfBbRr2YRZd5/NLf06eU4GoD5pT2rrStN6gIhRU0I4HegMLDTGnF++xeWPwHdA31AEJ0caMXcEhSWVyxAXlhRW6revy8wfXweIXbkubp1xa6Wd2eqZeoy/anyNF+vqBnZ7t+vNjMEzWDN0DTMGz6j0GhVdQJ5+5o/sDbu5/OVvmPjTBu46ryvT7jqTbm2aqb/bV7VN7/R17EWfv2Oq7TI6eIA7EYwBNgP9rLUbQxGYJ+oygs5jO7N2z9ojvp/YMpE1Q9cc0UfvbReLr90/dekugiPHEILVBVST0jLLa3NXMnbOChKaN2b0oBT6dTm2PMAq/d0VFzf1d9espq40Xz5Tff5B4XeXkTGmpTHmn7innA4EpgIfG2POD1yY4qvapk7WeeaPjwPEdZ3t5E0XkDctnLquf9iwM59B//yOFz9bzuWnHMfHw845lAzgiAFSVxtL+pV5dJ7YLyilPMKhVEhA1NSV5st6AO2w5qiaZhn9DLwG/NlaWwJ8aoxJAV4zxqyz1t4YkgijgCvXxYi5I3DluEhq657+WNe74YwBGcxZPeeIO+yKfvs6X6h9nL3kz2ynmvZw9qZ8Rl1KbFhrmb5wEyPfX4wBxg5K4apT23sI4FB/tysB+t8B+Q3A1itg3bKZAS3lES6lQmrkcsE998CPP7qf9+0LL73k+926tyU5/B1vCOWmQVGopjGEc6y1L5QnAwCstVnW2v7AF8EPzXnBWIXr78ri2u6wfV1cVqGmOfqePoeMARk0adCk0ms0adDE7y06vWnh+NoKmr8ui5RnX+G+ydnYBhsYfVNrz8kAKvV3jxhQkQzw6n2Cca6OcrmgXz/46qtDBe6+/NL9vWD16/uz1qMi3g8+cM96+uAD9/OpUzUm4aVqE0JNYwXW2jeCE074CNYqXF8rd3pKSDUNstZ1L+XqEg3g8XNYvmO5T5+Dt7xp4fjSCnon80cG/SOb3bs7savBv1la+keunnpO9dVP77/l4ACpK+FQMqjtfYJ1ro4aMcKdBKoqKAheF44/9YeGDj0y3oICGDSo5kVzcpDKX1cjUHdvdf2fvq4JyZ9pmp4STXWfw90f313rbKe68KaF480xB0pKefqjpfx16lbKKCSn8f3sbTgFa0rd/47PXILrtPb0//tplT/jebfi+ng8pKWRtC8OU6X3osbWlo+zY+ramguZmuIP1pRRf+oP/fCD5++XldU+JqGZTYDKX1crUHdvde1rrykh1TSLB2ruo/dVdZ/D9rztQbm7rW2MxJtjVuTuY+jELJZs2Qtx37LFjsGaA5XjtDmM6A755lAdloOf8bS7mZHdhIy+fZjT4EfySwurjeXQB+X7alxvztVRSUmeF5yB+2IdzPcN5oK1qmMSWkl9kKMtBGPMQGPMMmPMSmPMw07GUlVd796qdkHcknRLnbpwwqU7obrPoVV8q6Dc3XrTwqnumFPanMK/56/lile+IWdvIW/8NpXk7kvAFFWOswySc6i+S8jmwNq1JE3+mvlvQtpx59Xe2vI0OyYvDy65pNq7zWAvuvNbRoa7+qkn3pTIDpbq7ub79PHu96uOSWhm00G1rkMI2hsbUx9YDlwEbAR+Am601i6p7ndCuQ6hLvPlq/ud8VeN5x3XO17VB6pQl3n+gZzN5M053TrjVkfXExxu675CHpzq4stl2zive2ueuy6JNs2bHBl/GTQtgflvugeNZ55YOSmYMkhbBjMmVXzDyxIL1ZVvAHdRuEi922zfHjZvPvL7FWUpQq2mdQrgHkQ+fByhcWN3UissrH5dQ22lN6JAIEpXBFsfYKW1drW1tgiYCKQ7GE8ldbl7q66b5x3XOz6vtPV1cDhY+yRU9zlcd/J1YXN3+/mSXC4dO4/vVu1gVHpP3r61N22aN/Ec/962zH8TknIhY647OZgy9+tUJItRXx724t5OefQ0O6ZC1bvNSOqvTk0NrwqvNd3NJyW5q7Gmp7sv5unp7umy339feUxi/Hj361R8/p06hdc5OsmbjZeD8QCuA/512PNbgFc9HDcEyAQyO3bsWOdNpkMhcUyi5QmOeCSOTazT62XnZNv0Cek2cWyiTZ+QXuNm8WkT0g5uSl/xME8Ymz7Btw3Fg71BfSDlHSi2j0xz2U4PzbKXjv3aLs/ZW/svVdlEPjsBm35TfZv4XDubPqytzU44bEN2XzZlr3hd8PxITPT4/o5uWO+NcIs3MbHmz7c2ns4nLs79CJdzDAIg03pxXXZyUNnT7dQR/VfW2nHAOHB3GQU7KH8EujS1L4PDgRhziIiFUuVcG3czbGIWa3bk8YdzunDfxSfSuEH92n+xSlnnpORkZjxYXmbB5YI3+oOpQ3nrite95BLIyan8s8PvNmsrFR1u/N27ORjx+LNPsqfPv7AQzjsPWrQIj3N0kJNdRhuBDoc9Px53vaSIVdc1AIEQiCmMYb9QCncdor/PXck1r82noLiUd+/syyOX9fAuGVSorsyCv1suVrexzOFJJRSVPwPdJRVOFV793Se5us+/4tzC4Rwd5GRC+Ak4wRjT2RjTCBgMfOBgPH4L9KwRX1ZKByIZhcvMpups3JXPjeO+5/lPlnFJr7bMHnoO/bu2Cuyb+Hvxqy2pBHvXNW92L4tkgUjakTReEOrxJm/6lYL1AC7DPdNoFfBYbceffvrpAe9bC1fZOdk2/qn4g+MC5glj45+Kr7FP35cxB08CNQ4RDNN/3mh7jZxte46cbadmbrBlZWVOh1Q3we6TT0s79Nq+joOEm+xs9/kkJrq/BuIzCrcxkZoEMFa8HENwNCH4+oilhODExbkuSSjYducX2bv/97Pt9NAse+1r39r1O/IciyVgsrPdF+jERPfXQF6M/B10DRe1XQz9SRbB/PwDKYDJ3duEoJXKYcqJ7htf91QOtu9X7+C+SVnk7jvA8ItO5K7zutKgfhRUWwnmSlx/B13DRU2D76NG+beyONgroQPFgZ3mouD/rujkVJ2bQO5O5ssYyOHHXvm/qxk+9WtufON7GjWox3t39efuC07wmAyiZj+BQPF30NUpVfvKMzOrvxgGc2VxOK0RcWC8w7GVynURSzumhcPOYv7wJf7Dj61f1p7WRffTyHbj4l7xjLn+LOIbe27IRvpnFDQ17V4WjjytPq5Xr3JROji0ajw7Ozgri8Ntt7YAxhMJK5WlBmFf56YWvkxhHTF3BPlF+cSXXMpxB8ZS37ZmW6On2NpgbLXJwNf3iCj+3qWG0zRRb3i64y8tdScFTy2dYN05h1tNI39nVNWBxhDCWCCrloaaL2Mg2ZvW0KpoBE3L+lBQbwE7Go2l1OwiOzfx0Ot5qNMU7tNk6yRaK2/WtJOZp75ygIQE6N37yJZORob7M6l65+xvt5gDffYN5KSDAAAPEElEQVS1CvF4hxKCBIW3q7bnLM2lwa6RxJXVZ2fDf7Kv/iwwttKx1a2g7t2ud0BXhoeFSFvJ7I3aklx1A+G9e3s+Z0+rp2++2f+tM6NlQN4PGkOQoKitf7+gqJSnPlrCO9+vJ7FVQxYW3MO+suUej62u8uu5nc7lp80/RdcYQnWVN+Pi3HfMkbhPcHq6e4Gcp/GAGTP87ysPVF97uI0hBJDGEMRRNY2BLNq4h8tfmcc736/nzrM688mwC5g3ZHK14yXVdQ2t3bM2osdZPKquampBQeSuPK6tK8bfvvJA9f070GcfbtRCkJApLbP88+tVjP50Oa2aNebFG5I5s1vtpSfqsjdExKp6l+qJt3s0hIvaWgj+ioH9DPylFoKElU27C7jxje95bvYyLu6ZwOxhZ3uVDMDZooEhV/UuNS7uyGOcHuj0VbDXRkRafaIwpoQgQfd+1iYGjv2axZv28ML1yfz9N6fRsmkjr38/0qfg+uzwaaMXXRT5F7tgd8VE6mK8MKQuIwmaPQXFjHz/F97P2sxpHVsydtCpdDy2qdNhRZYoHugMqEhbjBdi3nYZadppFAjGXsr++mH1Du6bnE3O3kLuvfBE/jwgSuoQhVq4bVATriKlPlGYUwshwoVb+YaikjLGfr6c179aRcdjmjJmUAqndTw65HGIyCEaVI4R4VS+YdW2/Vz7+nxe+3IVN5zegQ/vOVvJQCSCqMsowoVD+QZrLe/+sJ6/fbiEJg3r84+bT2Ngr+NC9v4iEhhKCBHO2xIRwbJ9/wEefs/F50u3cvYJrXjh+mQSWjQJyXuLSGApIUS4jAEZzFk954gxhFDM0Z/761YemJrN3sISRlxxMrf1T6RePQ+rbEUkIjgyhmCMed4Y86sxxmWMmW6MaelEHNHAiTn6hcWljHz/F24b/xPHxjfmg7+cyR1ndVYyEIlwjswyMsZcDHxhrS0xxjwLYK19qLbf0ywj5y3evIehE7NYuXU/d5zVmQcu6U6ThvWdDktEahDW6xCstZ8e9vR74Don4hDvlZVZ3pi3mhc+XcbRTRvx3zv6cPYJrZ0OS0QCKBzGEG4HJlX3Q2PMEGAIQMeOHUMVkxxm8+4Chk/O5rvVOxjYsy1PX3MKR8d7X3pCRCJD0BKCMeZzoK2HHz1mrX2//JjHgBLg3epex1o7DhgH7i6jIIQqNZiZvZnHpi+ipMzy3LVJXJ96PMZTeWYRiXhBSwjW2gtr+rkx5nfAFcAFNpKWS8eIfYXFPP7+YqYt3MSpHVsydlAKnY6NdzosEQkiR7qMjDEDgYeAc621+U7EINXLXLuTYZOy2LKnkGEXnsBfBnRTHSKRGODUGMKrQGPgs/Luh++ttX90KBYpV1xaxstzVvD3uSs5/uimTP7DGZzeSaUnRGKFU7OMujnxvlK9NdvzGDYpi+wNu7n+9ON5PK0nzRqHw5wDEQkV/R8f46y1TPxpA6NmLqFRg3q8ftNpXHqK6hCJxCIlhBi2M6+Ih99z8emSXM7sdiwvXp9C26NUh0gkVikhxKivlm/j/inZ7Mkv5q+X9+D2M1V6QiTWKSHEmMLiUp75+FfGz19L94Tm/Of2PvQ4roXTYYlIGFBCiCFLNu9l2KSFLM/dz21nJvLQwJNUh0hEDlJCiAFlZZY3v1nD858s46imDfnP7X0450TVIRKRypQQotyWPQXcPyWbb1fu4OKTE3jm2iSOUR0iEfFACSGKfbRoC49MW0RxaRnPXHMKg3p3UB0iEamWEkIU2n+ghCc+WMzUBRtJ7uCuQ9S5leoQiUjNlBCizIJ1O7l3UjYbd+Vzz/nduPuCE2ioOkQi4gUlhChRUlrGy1+s5NUvVtCuZRyT/3AGqYnHOB2WiEQQJYQosLa8DlHWht1cc1p7nkzrSfMmDZ0OS0QijBJCBLPWMjlzA0/OXEKDeoZXbjyVK5PbOR2WiEQoJYQItSuviIenufhkcS5ndDmWF29Ipl3LOKfDEpEIpoQQgb4ur0O0K7+IRy49id+f3UV1iETEb0oIEaSwuJTnZi/jrW/X0K1NM96+rTc92x3ldFgiEiWUECLE0i17GTYxi2W5+/jdGZ145LIeqkMkIgHlaEIwxtwPPA+0ttZudzKWcFVWZnnr2zU8N3sZLeIa8vZtvRnQvY3TYYlIFHIsIRhjOgAXAeudiiHc5e4tZPjkbL5ZuZ0LeyTw7LWncGyzxk6HJSJRyskWwhjgQeB9B2MIW7N/2cLD0xZRWFzKU1f34jd9OqoOkYgElSMJwRiTBmyy1mbrIlfZ/gMlPPnBYqYs2EjS8UcxZlAKXVs3czosEYkBQUsIxpjPgbYefvQY8ChwsZevMwQYAtCxY8eAxReOfl6/i3snZbF+Zz5/HtCVoRecSKMGqkMkIqFhrLWhfUNjTgHmAPnl3zoe2Az0sdbm1PS7qampNjMzM8gRhl5JaRmvfLGSV+eupG2LJowZlEKfzqpDJCKBYYxZYK1Nre24kHcZWWsXAQenyRhj1gKpsTrLaN0Odx2ihet3c/Wp7XkyvSctVIdIRBygdQgOsdYyZcFGnvxgMfXqGV6+8VTSVIdIRBzkeEKw1iY6HUOo7cor4tHpi/j4lxz6dj6G0YNSaK86RCLiMMcTQqz5ZsV2hk/JYmdeEQ8NPIkh53ShvuoQiUgYUEIIkcLiUp7/ZBlvfrOGrq3jefN3venVXnWIRCR8KCGEwLKcfQyduJBfc/ZxS79OPHpZD+IaqQ6RiIQXJYQgKiuzjJ+/lmdm/0qLJg1469ZUzj8pwemwREQ8UkIIkq17Cxk+JZt5K7Zz/kltePbaJFo3Vx0iEQlfSghBMPuXHB6Z5qKguJS/XdWLm/qqDpGIhD8lhADKO1DCqJlLmJS5gV7tWzB20Kl0a6M6RCISGZQQAiRrw26GTVzIup35/Om8rgy7UHWIRCSyKCH4qaS0jNe+XMVLc1bQtkUTJvy+H/26HOt0WCIiPlNC8MOGnfncOymLzHW7SE9px6j0XhwVpzpEIhKZlBDqwFrLtJ838fgHizHAS4NTSE9p73RYIiJ+UULw0Z78Yh6dsYgPXVvok3gMowclc/zRTZ0OS0TEb0oIPpi/cjvDp2Szbd8BHhzYnT+c01V1iEQkaigheOFASSkvfrqcN+atpnOreKb/6UxOOV51iEQkuigh1GJF7j6GTsxiyZa93NS3I49d3oOmjfSxiUj00ZWtGtZa/vPdOv7vo6U0a9yAf/02lQtPVh0iEYleSggebN1XyANTXHy1fBsDurfmueuSVYdIRKKeEkIVny3J5aH3XOQdKCEjvSc39+ukOkQiEhMcSwjGmLuBvwAlwIfW2gedigUgv6iEjFlLmfDjenq2a8FLg1Po1qa5kyGJiISUIwnBGDMASAeSrLUHjDFtnIijQvaG3QyblMXaHXn84dwuDL+ou+oQiUjMcaqFcBfwjLX2AIC1dqsTQZSWWV7/ciVjP19B6+aNeffOvvTv2sqJUEREHOdUQjgRONsY8xRQCNxvrf0plAFs2JnPfZOz+GntLq5Mbsff0ntxVFPVIRKR2BW0hGCM+Rxo6+FHj5W/79FAP6A3MNkY08Vaaz28zhBgCEDHjh39jstay4ysTYycsRiAMYOSuSqlvQaORSTmBS0hWGsvrO5nxpi7gGnlCeBHY0wZ0ArY5uF1xgHjAFJTU49IGL7Yk1/MX9//hZnZm+mdeDSjb0ihwzGqQyQiAs51Gc0Azge+NMacCDQCtgfzDb9btYPhk7PYuu8AD1zSnT+eqzpEIiKHcyohvAW8ZYz5BSgCfuepuyhQXpmzgtGfLyfx2Hjeu6s/yR1aBuutREQiliMJwVpbBNwcqvdLbBXP4N4dGXGF6hCJiFQnJq6OVya348rkdk6HISIS1rT6SkREACUEEREpp4QgIiKAEoKIiJRTQhAREUAJQUREyikhiIgIoIQgIiLlTBArRgScMWYbsM7pOOqoFUGu1xTmYvn8Y/ncIbbPP1zOvZO1tnVtB0VUQohkxphMa22q03E4JZbPP5bPHWL7/CPt3NVlJCIigBKCiIiUU0IInXFOB+CwWD7/WD53iO3zj6hz1xiCiIgAaiGIiEg5JYQQMsZcb4xZbIwpM8ZEzMwDfxhjBhpjlhljVhpjHnY6nlAyxrxljNlavjNgTDHGdDDGzDXGLC3/mx/qdEyhZIxpYoz50RiTXX7+TzodkzeUEELrF+Aa4GunAwkFY0x94O/ApcDJwI3GmJOdjSqkxgMDnQ7CISXAcGttD6Af8OcY+7c/AJxvrU0GUoCBxph+DsdUKyWEELLWLrXWLnM6jhDqA6y01q4u3zZ1IpDucEwhY639GtjpdBxOsNZusdb+XP7f+4ClQHtnowod67a//GnD8kfYD9gqIUgwtQc2HPZ8IzF0URA3Y0wicCrwg7ORhJYxpr4xJgvYCnxmrQ3784+JPZVDyRjzOdDWw48es9a+H+p4HGY8fC/s75IkcIwxzYD3gGHW2r1OxxNK1tpSIMUY0xKYbozpZa0N6/EkJYQAs9Ze6HQMYWQj0OGw58cDmx2KRULMGNMQdzJ411o7zel4nGKt3W2M+RL3eFJYJwR1GUkw/QScYIzpbIxpBAwGPnA4JgkBY4wB3gSWWmtHOx1PqBljWpe3DDDGxAEXAr86G1XtlBBCyBhztTFmI3AG8KEx5hOnYwoma20J8BfgE9yDipOttYudjSp0jDETgO+A7saYjcaYO5yOKYTOBG4BzjfGZJU/LnM6qBA6DphrjHHhvjH6zFo7y+GYaqWVyiIiAqiFICIi5ZQQREQEUEIQEZFySggiIgIoIYiISDklBBEflFfxXGOMOab8+dHlzzsZY2YbY3YbY8J+eqGIJ0oIIj6w1m4AXgeeKf/WM8A4a+064Hncc+9FIpISgojvxgD9jDHDgLOAFwGstXOAfU4GJuIP1TIS8ZG1ttgY8wAwG7i4vLS3SMRTC0Gkbi4FtgC9nA5EJFCUEER8ZIxJAS7CvRPYvcaY4xwOSSQglBBEfFBexfN13PX91+MeSH7B2ahEAkMJQcQ3vwfWW2s/K3/+GnCSMeZcY8w8YApwQXl100sci1KkDlTtVEREALUQRESknBKCiIgASggiIlJOCUFERAAlBBERKaeEICIigBKCiIiUU0IQEREA/h9N1aJaGbB7qQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1b430a37550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### 绘制逻辑回归决策边界\n",
    "def plot_decision boundary(X_train, y_train, params):\n",
    "    '''\n",
    "    输入：\n",
    "    X_train: 训练集输入\n",
    "    y_train: 训练集标签\n",
    "    params：训练好的模型参数\n",
    "    输出：\n",
    "    决策边界图\n",
    "    '''\n",
    "    # 训练样本量\n",
    "    n = X_train.shape[0]\n",
    "    # 初始化类别坐标点列表\n",
    "    xcord1 = []\n",
    "    ycord1 = []\n",
    "    xcord2 = []\n",
    "    ycord2 = []\n",
    "    # 获取两类坐标点并存入列表\n",
    "    for i in range(n):\n",
    "        if y_train[i] == 1:\n",
    "            xcord1.append(X_train[i][0])\n",
    "            ycord1.append(X_train[i][1])\n",
    "        else:\n",
    "            xcord2.append(X_train[i][0])\n",
    "            ycord2.append(X_train[i][1])\n",
    "    # 创建绘图\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    # 绘制两类散点，以不同颜色表示\n",
    "    ax.scatter(xcord1, ycord1,s=32, c='red')\n",
    "    ax.scatter(xcord2, ycord2, s=32, c='green')\n",
    "    # 取值范围\n",
    "    x = np.arange(-1.5, 3, 0.1)\n",
    "    # 决策边界公式\n",
    "    y = (-params['b'] - params['W'][0] * x) / params['W'][1]\n",
    "    # 绘图\n",
    "    ax.plot(x, y)\n",
    "    plt.xlabel('X1')\n",
    "    plt.ylabel('X2')\n",
    "    plt.show()\n",
    "    \n",
    "plot_logistic(X_train, y_train, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.datasets.samples_generator import make_classification\n",
    "\n",
    "class logistic_regression():\n",
    "    def __init__(self):\n",
    "        pass\n",
    "\n",
    "    def sigmoid(self, x):\n",
    "        z = 1 / (1 + np.exp(-x))\n",
    "        return z\n",
    "\n",
    "    def initialize_params(self, dims):\n",
    "        W = np.zeros((dims, 1))\n",
    "        b = 0\n",
    "        return W, b\n",
    "\n",
    "    def logistic(self, X, y, W, b):\n",
    "        num_train = X.shape[0]\n",
    "        num_feature = X.shape[1]\n",
    "\n",
    "        a = self.sigmoid(np.dot(X, W) + b)\n",
    "        cost = -1 / num_train * np.sum(y * np.log(a) + (1 - y) * np.log(1 - a))\n",
    "\n",
    "        dW = np.dot(X.T, (a - y)) / num_train\n",
    "        db = np.sum(a - y) / num_train\n",
    "        cost = np.squeeze(cost)\n",
    "\n",
    "        return a, cost, dW, db\n",
    "\n",
    "    def logistic_train(self, X, y, learning_rate, epochs):\n",
    "        W, b = self.initialize_params(X.shape[1])\n",
    "        cost_list = []\n",
    "\n",
    "        for i in range(epochs):\n",
    "            a, cost, dW, db = self.logistic(X, y, W, b)\n",
    "            W = W - learning_rate * dW\n",
    "            b = b - learning_rate * db\n",
    "\n",
    "            if i % 100 == 0:\n",
    "                cost_list.append(cost)\n",
    "            if i % 100 == 0:\n",
    "                print('epoch %d cost %f' % (i, cost))\n",
    "\n",
    "        params = {\n",
    "            'W': W,\n",
    "            'b': b\n",
    "        }\n",
    "        grads = {\n",
    "            'dW': dW,\n",
    "            'db': db\n",
    "        }\n",
    "\n",
    "        return cost_list, params, grads\n",
    "\n",
    "    def predict(self, X, params):\n",
    "        y_prediction = self.sigmoid(np.dot(X, params['W']) + params['b'])\n",
    "\n",
    "        for i in range(len(y_prediction)):\n",
    "            if y_prediction[i] > 0.5:\n",
    "                y_prediction[i] = 1\n",
    "            else:\n",
    "                y_prediction[i] = 0\n",
    "\n",
    "        return y_prediction\n",
    "\n",
    "    def accuracy(self, y_test, y_pred):\n",
    "        correct_count = 0\n",
    "        for i in range(len(y_test)):\n",
    "            for j in range(len(y_pred)):\n",
    "                if y_test[i] == y_pred[j] and i == j:\n",
    "                    correct_count += 1\n",
    "\n",
    "        accuracy_score = correct_count / len(y_test)\n",
    "        return accuracy_score\n",
    "\n",
    "    def create_data(self):\n",
    "        X, labels = make_classification(n_samples=100, n_features=2, n_redundant=0, n_informative=2,\n",
    "                                        random_state=1, n_clusters_per_class=2)\n",
    "        labels = labels.reshape((-1, 1))\n",
    "        offset = int(X.shape[0] * 0.9)\n",
    "        X_train, y_train = X[:offset], labels[:offset]\n",
    "        X_test, y_test = X[offset:], labels[offset:]\n",
    "        return X_train, y_train, X_test, y_test\n",
    "\n",
    "    def plot_logistic(self, X_train, y_train, params):\n",
    "        n = X_train.shape[0]\n",
    "        xcord1 = []\n",
    "        ycord1 = []\n",
    "        xcord2 = []\n",
    "        ycord2 = []\n",
    "        for i in range(n):\n",
    "            if y_train[i] == 1:\n",
    "                xcord1.append(X_train[i][0])\n",
    "                ycord1.append(X_train[i][1])\n",
    "            else:\n",
    "                xcord2.append(X_train[i][0])\n",
    "                ycord2.append(X_train[i][1])\n",
    "        fig = plt.figure()\n",
    "        ax = fig.add_subplot(111)\n",
    "        ax.scatter(xcord1, ycord1, s=32, c='red')\n",
    "        ax.scatter(xcord2, ycord2, s=32, c='green')\n",
    "        x = np.arange(-1.5, 3, 0.1)\n",
    "        y = (-params['b'] - params['W'][0] * x) / params['W'][1]\n",
    "        ax.plot(x, y)\n",
    "        plt.xlabel('X1')\n",
    "        plt.ylabel('X2')\n",
    "        plt.show()\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    model = logistic_regression()\n",
    "    X_train, y_train, X_test, y_test = model.create_data()\n",
    "    print(X_train.shape, y_train.shape, X_test.shape, y_test.shape)\n",
    "    cost_list, params, grads = model.logistic_train(X_train, y_train, 0.01, 1000)\n",
    "    print(params)\n",
    "    y_train_pred = model.predict(X_train, params)\n",
    "    accuracy_score_train = model.accuracy(y_train, y_train_pred)\n",
    "    print('train accuracy is:', accuracy_score_train)\n",
    "    y_test_pred = model.predict(X_test, params)\n",
    "    accuracy_score_test = model.accuracy(y_test, y_test_pred)\n",
    "    print('test accuracy is:', accuracy_score_test)\n",
    "    model.plot_logistic(X_train, y_train, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Installation\\anaconda\\install\\lib\\site-packages\\sklearn\\linear_model\\logistic.py:433: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.\n",
      "  FutureWarning)\n",
      "D:\\Installation\\anaconda\\install\\lib\\site-packages\\sklearn\\utils\\validation.py:761: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0, 1, 1, 0, 1, 1, 0, 0, 1, 0])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "clf = LogisticRegression(random_state=0).fit(X_train, y_train)\n",
    "y_pred = clf.predict(X_test)\n",
    "y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       1.00      1.00      1.00         5\n",
      "           1       1.00      1.00      1.00         5\n",
      "\n",
      "   micro avg       1.00      1.00      1.00        10\n",
      "   macro avg       1.00      1.00      1.00        10\n",
      "weighted avg       1.00      1.00      1.00        10\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(classification_report(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
